## Replication file for:
## ``Off-Cycle and Out of Office: Election Timing and the Incumbency Advantage"
## Justin de Benedictis-Kessner

rm(list=ls())
library(rdrobust)
library(rdd)
library(plyr); library(dplyr)
library(stringr); library(reshape2)
library(tools)
library(ggplot2)
library(xtable)
set.seed(02139)
setwd("Replication") # Set working directory to folder with all replication files
source("rd.export.R")


#######################
#### Read the Data ####
#######################
# Elections data from de Benedictis-Kessner & Warshaw, 2016, ``Mayoral Partisanship and Municipal Fiscal Policy", Journal of Politics 78(4): 1124-1138. Updated with additional cities under 75,000 in population from Ferreira & Gyourko and to fill in temporal holes in existing cities's data, and merged with ICMA institutions data.
data <- read.csv(file="mayoral_elecs_updated20170413.csv",stringsAsFactors=F) # 9132 elections

## Media data scraped from Library of Congress:
media <- read.csv("city_papers_20k.csv",stringsAsFactors = F)

#############################
#### Set up incumbent VS ####
#############################

###########
##### Construct incumbency and VS measures:
###########
# fix month var:
for(i in 1:nrow(data)){
  if(!is.na(data$elecdate[i])){
    data$month[i] <- gsub("\\d{4}-(\\d{2})-\\d{2}","\\1",x=data$elecdate[i])
  }
}
sum(is.na(data$month)) # 2681 NAs

data$use2 <- NA

data$mayor_next_VS <- NA
data$runnerup_next_VS <- NA
data$mayor_run_next <- NA
data$runnerup_run_next <- NA
data$elec_index_next <- NA
data$YearData_next <- NA
nextelection <- NULL
for(i in 1:nrow(data)){
  thiscity <- data[which(data$fips==data$fips[i]),]
  futureyears <- thiscity[which(thiscity$YearData > data$YearData[i]),]
  if(nrow(futureyears)>0){
    nextelection <- futureyears[order(futureyears$YearData,futureyears$month, decreasing=F)[1],]
    data$elec_index_next[i] <- nextelection$elec_index
    data$YearData_next[i] <- nextelection$YearData
    if(!is.na(data$mayor_name_final[i])){
      if(data$YearData_next[i]-data$YearData[i] <= 6){ # only use data for "next election" if reasonable time frame
      if(data$mayor_name_final[i] %in% nextelection[,c("mayor_name_final","runnerup_name_final")]){
        data$mayor_run_next[i] <- 1
        data$use2[i] <- 1
      } 
      if(!(data$mayor_name_final[i] %in% nextelection[,c("mayor_name_final","runnerup_name_final")])){
        data$mayor_run_next[i] <- 0
      }
    }
    if(!is.na(data$runnerup_name_final[i])){
      if(data$runnerup_name_final[i] %in% nextelection[,c("mayor_name_final","runnerup_name_final")]){
        data$runnerup_run_next[i] <- 1
        data$use2[i] <- 1
      } 
      if(!(data$runnerup_name_final[i] %in% nextelection[,c("mayor_name_final","runnerup_name_final")])){
        data$runnerup_run_next[i] <- 0
      }
    }
    }
    if(is.na(data$mayor_name_final[i])){
      data$mayor_run_next[i] <- NA
    }
    if(is.na(data$runnerup_name_final[i])){
      data$runnerup_run_next[i] <- NA
    }
    if(!(data$mayor_name_final[i] %in% nextelection[,c("mayor_name_final","runnerup_name_final")]) & !(data$runnerup_name_final[i] %in% nextelection[,c("mayor_name_final","runnerup_name_final")])){
      data$use2[i] <- 0
    }
    if(is.na(data$runnerup_name_final[i]) & is.na(data$runnerup_name_final[i])){
      data$use2[i] <- 0
    }
  }
  if(nrow(futureyears)<1 | (data$YearData_next[i]-data$YearData[i] > 6)){
    data$mayor_run_next[i] <- NA
    data$runnerup_run_next[i] <- NA
    data$use2[i] <- 0
  }
  
  if(data$use2[i]==1){
    thiscity <- data[which(data$fips==data$fips[i]),]
    futureyears <- thiscity[which(thiscity$YearData > data$YearData[i]),]
    
    if(nrow(futureyears)>0){
      nextelection <- futureyears[order(futureyears$YearData, decreasing=F)[1],]
      
      if(!is.na(nextelection$mayor_votes_final) & !is.na(nextelection$runnerup_votes_final)){
        if(!is.na(data$mayor_name_final[i])){
          if(!is.na(nextelection$mayor_name_final) & data$mayor_name_final[i] == nextelection$mayor_name_final){
            data$mayor_next_VS[i] <- nextelection$mayor_votes_final/(nextelection$mayor_votes_final+nextelection$runnerup_votes_final)
          }
          if(!is.na(nextelection$runnerup_name_final) & data$mayor_name_final[i] == nextelection$runnerup_name_final){
            data$mayor_next_VS[i] <- nextelection$runnerup_votes_final/(nextelection$mayor_votes_final+nextelection$runnerup_votes_final)
          }
        }
        if(is.na(data$mayor_name_final[i])){
          data$mayor_next_VS[i] <- NA
        }
        
        if(!is.na(data$runnerup_name_final[i])){
          if(!is.na(nextelection$mayor_name_final) & data$runnerup_name_final[i] == nextelection$mayor_name_final){
            data$runnerup_next_VS[i] <- nextelection$mayor_votes_final/(nextelection$mayor_votes_final+nextelection$runnerup_votes_final)
          }
          if(!is.na(nextelection$runnerup_name_final) & data$runnerup_name_final[i] == nextelection$runnerup_name_final){
            data$runnerup_next_VS[i] <- nextelection$runnerup_votes_final/(nextelection$mayor_votes_final+nextelection$runnerup_votes_final)
          } 
        }
        if(is.na(data$runnerup_name_final[i])){
          data$runnerup_next_VS[i] <- NA
        }
      }
      if(is.na(nextelection$mayor_votes_final) | is.na(nextelection$runnerup_votes_final)){
        data$mayor_next_VS[i] <- NA
        data$runnerup_next_VS[i] <- NA
      }
    }
  }
  if(data$use2[i]==0){
    data$mayor_prev_VS[i] <- NA
    data$runnerup_prev_VS[i] <- NA
  }
}

# summary(data$mayor_next_VS);sum(!is.na(data$mayor_next_VS)) # 65.4%, 4203 with data (out of 9186)!
# summary(data$runnerup_next_VS) # 35.3%
# summary(data$mayor_run_next) # 62.16%
# summary(data$runnerup_run_next) # 19.01%

data$mayor_voteshare <- data$mayor_votes_final/(data$mayor_votes_final + data$runnerup_votes_final)
data$runnerup_voteshare <- data$runnerup_votes_final/(data$mayor_votes_final + data$runnerup_votes_final)
data$turnout <- data$mayor_votes_final + data$runnerup_votes_final

# checking elec_index matching:
# head(table(data$elec_index_next)[order(table(data$elec_index_next),decreasing = T)],10)
# fix elections with multiple elections in years: 965 1708 3270 4790 5800 
data$elec_index_next[which(data$elec_index==4784)] <- 4785 # fixing from 4790
data$elec_index_next[which(data$elec_index==5797)] <- 5796 # fixing from 5800
data$elec_index_next[which(data$elec_index==3268)] <- 3269 # fixing from 3270
data$elec_index_next[which(data$elec_index==964)] <- 9398 # fixing from 965
data$elec_index_next[which(data$elec_index==1705)] <- 1704 # fixing from 1708



## election timing:

natyears <- c(seq(from=1940,to=2014,by=2))
presyears <-c(seq(from=1940,to=2012,by=4))
midyears <- c(seq(from=1942,to=2014,by=4))

data$concurrent <- NA
data$concurrent_pres <- NA
data$concurrent_mid <- NA
for(i in 1:nrow(data)){
  if(!is.na(data$YearData[i]) & !is.na(data$month[i])){
  	
    if(data$YearData[i] %in% natyears & data$month[i]==11){
      data$concurrent[i] <- 1
    }
    if(data$YearData[i] %in% natyears & data$month[i]!=11){
      data$concurrent[i] <- 0
    }
    if(!(data$YearData[i] %in% natyears)){
      data$concurrent[i] <- 0
    }
    
    if(data$YearData[i] %in% presyears & data$month[i]==11){
      data$concurrent_pres[i] <- 1
    }
    if(data$YearData[i] %in% presyears & data$month[i]!=11){
      data$concurrent_pres[i] <- 0
    }
    if(!(data$YearData[i] %in% presyears)){
      data$concurrent_pres[i] <- 0
    }
    
    if(data$YearData[i] %in% midyears & data$month[i]==11){
      data$concurrent_mid[i] <- 1
    }
    if(data$YearData[i] %in% midyears & data$month[i]!=11){
      data$concurrent_mid[i] <- 0
    }
    if(!(data$YearData[i] %in% midyears)){
      data$concurrent_mid[i] <- 0
    }
    
  }
  if(!is.na(data$YearData[i]) & is.na(data$month[i])){ # if only year of data:
  	if(!(data$YearData[i] %in% natyears)){
      data$concurrent[i] <- 0
    }
    if(!(data$YearData[i] %in% presyears)){
      data$concurrent_pres[i] <- 0
    }
    if(!(data$YearData[i] %in% midyears)){
      data$concurrent_mid[i] <- 0
    }
  }
  if(is.na(data$month[i]) & is.na(data$elecdate[i]) & is.na(data$YearData[i])){
    data$concurrent[i] <- NA
    data$concurrent_pres[i] <- NA
    data$concurrent_mid[i] <- NA
  }
}
# summary(data$concurrent);table(data$concurrent) # 7.39%, 615 elections
# summary(data$concurrent_pres);table(data$concurrent_pres) # 3.28%, 288 elections

# create lead var
data$concurrent_next <- data$concurrent[match(data$elec_index_next,data$elec_index)]
# summary(data$concurrent_next); sum(data$concurrent_next==1,na.rm=T) # only 7.6% concurrent, or 544 elections! (2029 NAs)
data$concurrent_pres_next <- data$concurrent_pres[match(data$elec_index_next,data$elec_index)]
data$concurrent_mid_next <- data$concurrent_mid[match(data$elec_index_next,data$elec_index)]
# summary(data$concurrent_pres_next); sum(data$concurrent_pres_next==1,na.rm=T) # 3.4%, or 266 elections

## Construct measure of demshare in next election:
data$demshare_next <- data$demshare[match(data$elec_index_next,data$elec_index)]

for(i in 1:nrow(data)){
  if(!is.na(data$elec_index_next[i])){
    data$demshare_next[i] <- ifelse((data$mayor_party_final[which(data$elec_index==data$elec_index_next[i])]!="D" & data$runnerup_party_final[which(data$elec_index==data$elec_index_next[i])]!="D"),
                                    0, # change to zero if neither candidate was a Dem
                                    data$demshare_next[i])
  }
}

##########
### Manipulate into a candidate-level dataframe:
##########

data2 <- data[,c("fips","Name","STATEAB","JURIS","YearData","YearData_next","elec_index","elec_index_next","elecdate","month","population_est","mayor_name_final","runnerup_name_final","mayor_party_final","runnerup_party_final","mayor_voteshare","runnerup_voteshare","mayor_next_VS","runnerup_next_VS","mayor_run_next","runnerup_run_next")]
data3 <- melt(data2,id.vars=c("fips","Name","STATEAB","JURIS","YearData","YearData_next","elec_index","elec_index_next","elecdate","month","population_est","mayor_name_final","mayor_party_final"),measure.vars= c("mayor_voteshare","mayor_next_VS","mayor_run_next"))
data4 <- dcast(data3,fips + Name + STATEAB + JURIS + YearData + YearData_next + elec_index + elec_index_next + elecdate + month + population_est + mayor_name_final + mayor_party_final ~ variable)
names(data4)[grep("mayor_name_final",x=names(data4))] <- "cand_name"
names(data4)[grep("mayor_party_final",x=names(data4))] <- "cand_party"
names(data4)[grep("mayor_voteshare",x=names(data4))] <- "voteshare"
names(data4)[grep("mayor_next_VS",x=names(data4))] <- "next_VS"
names(data4)[grep("mayor_run_next",x=names(data4))] <- "run_next"


data3b <- melt(data2,id.vars=c("fips","Name","STATEAB","JURIS","YearData","YearData_next","elec_index","elec_index_next","elecdate","month","population_est","runnerup_name_final","runnerup_party_final"),measure.vars= c("runnerup_voteshare","runnerup_next_VS","runnerup_run_next"))
data4b <- dcast(data3b,fips + Name + STATEAB + JURIS + YearData + YearData_next + elec_index + elec_index_next + elecdate + month + population_est + runnerup_name_final + runnerup_party_final ~ variable)
names(data4b)[grep("runnerup_name_final",x=names(data4b))] <- "cand_name"
names(data4b)[grep("runnerup_party_final",x=names(data4b))] <- "cand_party"
names(data4b)[grep("runnerup_voteshare",x=names(data4b))] <- "voteshare"
names(data4b)[grep("runnerup_next_VS",x=names(data4b))] <- "next_VS"
names(data4b)[grep("runnerup_run_next",x=names(data4b))] <- "run_next"
data5 <- rbind(data4,data4b) # all set

# join in turnout, demshare
data5$turnout <- data$turnout[match(data5$elec_index,data$elec_index)]

# join back in institutions data:
data5$fog <- data$fog[match(data5$elec_index_next,data$elec_index)]
data5$partisan <- data$partisan[match(data5$elec_index_next,data$elec_index)]
data5$initiative <- data$initiative[match(data5$elec_index_next,data$elec_index)]
data5$referendum <- data$referendum[match(data5$elec_index_next,data$elec_index)]

data5$concurrent_next <- data$concurrent_next[match(data5$elec_index,data$elec_index)]
data5$concurrent_pres_next <- data$concurrent_pres_next[match(data5$elec_index,data$elec_index)]
data5$concurrent_mid_next <- data$concurrent_mid_next[match(data5$elec_index,data$elec_index)]

# summary(data5$concurrent_next) # still 7.6%, double the NAs from election-level data (3446)
# summary(data5$concurrent_pres_next) #  3.4%
# summary(data5$concurrent_mid_next) #  3.9%

# construct run/win combo dummy variable:
data5$run_and_win_next <- NA
for(i in 1:nrow(data5)){
  
  if(!is.na(data5$run_next[i]) & !is.na(data5$next_VS[i])){
    data5$run_and_win_next[i] <- ifelse(data5$run_next[i]==1 & data5$next_VS[i]>=0.5,1,0)  
  }
}
summary(data5$run_and_win_next) # 69.3% of candidates
# table(data5$run_and_win_next,data5$run_next) # shows that people who don't run are coded NA in run_and_win, need to change to 0
for(i in 1:nrow(data5)){
  if(!is.na(data5$run_next[i]) & data5$run_next[i]==0){
    data5$run_and_win_next[i] <- 0
  }
}
summary(data5$run_and_win_next) # 24.7% of candidates

#####
### Prep variables to join with media market data:
#####
data5$city <- gsub( " TOWNSHIP", "",data5$Name)
data5$city <- gsub( " TOWN", "",data5$city)
data5$city <- gsub( " CITY AND BOROUGH", "",data5$city)
data5$city <- gsub( " CITY AND COUNTY", "",data5$city)
data5$city <- gsub( " CITY-PARISH", "",data5$city)
data5$city <- gsub( " CITY", "",data5$city)
data5$city <- gsub( " VILLAGE", "",data5$city)
data5$city <- gsub( "ST ", "St. ",data5$city)
data5$city <- gsub( "FT ", "Fort ",data5$city)
data5$city <- gsub( "MT ", "Mount ",data5$city)
data5$city <- gsub( "WASHINGTON DC", "Washington",data5$city)
data5$city <- gsub( "WILKES BARRE", "WILKES-BARRE",data5$city)
data5$city <- gsub( "WINSTON SALEM", "WINSTON-SALEM",data5$city)
data5$city <- gsub( " MUNICIPALITY", "",data5$city)
data5$city <- tolower(data5$city)

data5$city <- gsub( "lexington-fayette urban co govt", "lexington",data5$city)
data5$city <- gsub( "lexington-fayette urban county", "lexington",data5$city)
data5$city <- gsub( "lexington-fayette urban county g", "lexington",data5$city)
data5$city <- gsub( "lexington-fayette urban county government",
                    "lexington",data5$city)

#####
### Merge with media data: 
#####
media$city <- tolower(media$city)
data6 <- merge(data5,media,by.x=c("city","STATEAB","YearData_next"),by.y=c("city","abb","years"),all.x=T)
data6 <- subset(data6,YearData>=1950)
data6$voteshare_adj <- data6$voteshare-0.5
data6$fog2 <- recode(data6$fog,"1=1;2=0;else=NA")

## calculate dailypaper difference variable (if it gained/lost paper):
data6$dailypaper_delta <- NULL
for(i in 1:nrow(data6)){
  thiscity <- data6[which(data6$fips==data6$fips[i]),]
  prevyears <- thiscity[which(thiscity$YearData<data6$YearData[i]),]
  if(nrow(prevyears)>0){
    prevyears <- prevyears[order(prevyears$YearData,decreasing=T),]
    data6$dailypaper_delta[i] <- data6$dailypaper[i]- prevyears$dailypaper[1]
  }
  if(nrow(prevyears)==0){
    data6$dailypaper_delta[i] <- NA
  }
}
table(data6$dailypaper_delta) 

data6$dailypaper_delta2 <- NULL
for(i in 1:nrow(data6)){
  thiscity <- data6[which(data6$fips==data6$fips[i]),]
  prevyears <- thiscity[which(thiscity$YearData<data6$YearData[i]),]
  if(nrow(prevyears)>0){
    prevyears <- prevyears[order(prevyears$YearData,decreasing=T),]
    data6$dailypaper_delta2[i] <- prevyears$dailypaper_delta[1]
  }
  if(nrow(prevyears)==0){
    data6$dailypaper_delta2[i] <- NA
  }
}
table(data6$dailypaper_delta2) # 68,22, 9140

# delta anticipatory variables:
data6$dailypaper_delta_lead <- NULL
data6$dailypaper_delta_lead2 <- NULL
for(i in 1:nrow(data6)){
  thiscity <- data6[which(data6$fips==data6$fips[i]),]
  futureyears <- thiscity[which(thiscity$YearData>data6$YearData[i]),]
  if(nrow(futureyears)>0){
    futureyears <- futureyears[order(futureyears$YearData,decreasing=F),]
    data6$dailypaper_delta_lead[i] <- futureyears$dailypaper_delta[1]
    if(length(unique(futureyears$YearData))>1){
      data6$dailypaper_delta_lead2[i] <- futureyears$dailypaper_delta[which(futureyears$YearData>futureyears$YearData[1])][1]
    }
    if(length(unique(futureyears$YearData))<2){
      data6$dailypaper_delta_lead2[i] <- NA
    }
  }
  if(nrow(futureyears)==0){
    data6$dailypaper_delta_lead[i] <- NA
    data6$dailypaper_delta_lead2[i] <- NA
  }
}


## election timing difference var:
data6$concurrent <- data$concurrent[match(data6$elec_index,data$elec_index)]
data6$concurrent_delta <- data6$concurrent_next - data6$concurrent


# delta anticipatory variables (in advance of a change in election timing):
data6$concurrent_delta_lead <- NULL
data6$concurrent_delta_lead2 <- NULL
for(i in 1:nrow(data6)){
  thiscity <- data6[which(data6$fips==data6$fips[i]),]
  futureyears <- thiscity[which(thiscity$YearData>data6$YearData[i]),]
  if(nrow(futureyears)>0){
    futureyears <- futureyears[order(futureyears$YearData,decreasing=F),]
    data6$concurrent_delta_lead[i] <- futureyears$concurrent_delta[1]
    if(length(unique(futureyears$YearData))>1){
      data6$concurrent_delta_lead2[i] <- futureyears$concurrent_delta[which(futureyears$YearData>futureyears$YearData[1])][1]
    }
    if(length(unique(futureyears$YearData))<2){
      data6$concurrent_delta_lead2[i] <- NA
    }
  }
  if(nrow(futureyears)==0){
    data6$concurrent_delta_lead[i] <- NA
    data6$concurrent_delta_lead2[i] <- NA
  }
}
table(data6$concurrent_delta_lead) # 282 about to go off-cycle, 196 about to go on-cycle
table(data6$concurrent_delta_lead2) # 266, 182

data6$concurrent_delta2 <- NULL
for(i in 1:nrow(data6)){
  thiscity <- data6[which(data6$fips==data6$fips[i]),]
  prevyears <- thiscity[which(thiscity$YearData<data6$YearData[i]),]
  if(nrow(prevyears)>0){
    prevyears <- prevyears[order(prevyears$YearData,decreasing=T),]
    data6$concurrent_delta2[i] <- prevyears$concurrent_delta[1]
  }
  if(nrow(prevyears)==0){
    data6$concurrent_delta2[i] <- NA
  }
}
sum(data6$concurrent_delta2==1,na.rm=T) 
sum(!is.na(data6$concurrent_delta2)) 

nrow(data) # 9131 elections
data6 <- data6[which(data6$cand_name != ""),] # 16120 obs
length(unique(data6$elec_index)) # 8974
nrow(data[!duplicated(data[,c("fips","YearData","mayor_votes_final")]),]) # 9128 unique election-years (multi elections in some years)
length(unique(data$fips)) # 1024 cities
length(unique(data6$fips)) # 1024 cities
nrow(unique(data6[which(data6$cand_name!=""),c("fips","cand_name")])) # 9919 unique candidates
nrow(unique(data6[which(data6$cand_name!=""),c("elec_index","cand_name")])) # 16086 individual obs

##################
#### Mapping: ####
##################

library(maptools) # lots of good tools for manipulating maps
library(rgdal) # better for opening shapefiles
library(RgoogleMaps) # googlemaps API tool
library(classInt) # color scales
library(RColorBrewer) # color scales
library(ggplot2); library(ggmap) # ggplot addition for fancy maps
library(maps);library(mapdata)

a <- data6[!duplicated(data6[,c("fips","city","STATEAB","concurrent")]),c("fips","city","STATEAB","concurrent","population_est","YearData")]
concurrenttab <- ddply(a,~fips,summarize,
                       max_timing = max(concurrent,na.rm=T))
a$concurrentmax <- concurrenttab$max_timing[match(a$fips,concurrenttab$fips)]
a <- a[!duplicated(a$fips),] # down to 1016 cities
a$fulladdress <- paste(a$city,a$STATEAB,sep=", ")
a$lon <- NA
a$lat <- NA
a[,c("lon","lat")] <- geocode(a$fulladdress,source = "google")




a$plotname <- NULL
for(i in 1:nrow(a)){
  a$plotname[i] <- ifelse(a$population_est[i]>=max(a$population_est[which(a$STATEAB==a$STATEAB[i])],na.rm=T),1,0)
  a$plotname[i] <- ifelse(a$population_est[i]>=500000,1,a$plotname[i])
}


simpleCap <- function(x) {
  s <- strsplit(x, " ")[[1]]
  paste(toupper(substring(s, 1,1)), substring(s, 2),
        sep="", collapse=" ")
}



# Problems with Alaska, Hawaii:
states <- readOGR("states_21basic","states")
states2 <- spTransform(states, CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))
states2@data$id <- rownames(states2@data)

alaska <- states2[which(states2$STATE_ABBR=="AK"),]
alaska <- elide(alaska, rotate=-50)
alaska <- elide(alaska, bb=bbox(alaska), scale=max(apply(bbox(alaska), 1, diff)) / 2.3)
alaska <- elide(alaska, shift=c(-2100000, -2500000))
proj4string(alaska) <- proj4string(states2)

# extract, then rotate & shift hawaii
hawaii <- states2[which(states2$STATE_ABBR=="HI"),]
hawaii <- elide(hawaii, rotate=-35)
hawaii <- elide(hawaii, shift=c(5400000, -1400000))
proj4string(hawaii) <- proj4string(states2)

states2 <- states2[!states2$STATE_ABBR %in% c("AK","HI"),]
states2 <- rbind(states2, alaska, hawaii)


# do the same for city points:
b <- data.frame(lon=a$lon,lat=a$lat,fips=a$fips,city=a$city,state=a$STATEAB,population=a$population_est,plotname=a$plotname,concurrent=a$concurrent,concurrentmax=a$concurrentmax,YearData=a$YearData)
b <- b[!is.na(b$fips),]

a <- data6[!duplicated(data6[,c("fips","city","STATEAB","concurrent")]),c("fips","city","STATEAB","concurrent","population_est","YearData")]
concurrenttab <- ddply(a,~fips,summarize,
                       min_timing = min(concurrent,na.rm=T))
b$concurrentmin <- concurrenttab$min_timing[match(b$fips,concurrenttab$fips)]

coordinates(b) <- ~ lon + lat
proj4string(b) <- CRS("+init=epsg:4326")
b <- spTransform(b,CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))

# adjusting anchorage + fairbanks coordinates:
ak_cities <- data.frame(b[which(b$state=="AK"),])
# locator(1)
ak_cities$lon <- c(-1404540,-1418058) # anchorage, fairbanks
ak_cities$lat=c(-2035876,-1879628) # anchorage, fairbanks
coordinates(ak_cities) <- ~ lon + lat
proj4string(ak_cities) <- proj4string(b)

b <- b[which(b$state!="AK"),]
b <- rbind(b,ak_cities[,!(names(ak_cities)=="optional")])


#####################
##### Figure 1: #####
#####################
pdf("cities_map1.pdf",height=10,width=13)
par(mar=c(0,0,0,0))
plot(states2,lwd=0.5,col="lightgrey",lty=1, border="white")
points(b,col="black",pch=".",cex=3)
text(b$lon[which(b$plotname==1)],b$lat[which(b$plotname==1)],labels=sapply(as.character(b$city[which(b$plotname==1)]),simpleCap),cex=0.8,pos=3,offset=0.2,col="black")
dev.off()
#####################



#####################
##### Figure 2: #####
#####################
pdf("elections_timeline.pdf",height=3,width=9)
ggplot(data) +
  geom_histogram(aes(x=YearData),binwidth=1) + 
  scale_x_continuous(breaks=seq(1950,2010,by=5),limits=c(1949,2015)) +
  theme_bw() +
  labs(x = "Year",
       y = "Count of elections") + 
  theme(axis.text = element_text(size=15),axis.title = element_text(size=15))
dev.off()
#####################


#### Placebo tests
data6$voteshare_last <- NA
data6$run_last <- NA
data6$run_and_win_last <- NA
for(i in 1:nrow(data6)){
  prevelec <- data6[which(data6$elec_index_next==data6$elec_index[i]),]
  if(nrow(prevelec)>0){ # if there was data for previous election
    if(!is.na(data6$cand_name[i]) & data6$cand_name[i] %in% prevelec$cand_name){ # if cand ran in last elec
      match <- prevelec[which(prevelec$cand_name==data6$cand_name[i]),]
      data6$run_last[i] <- 1
      data6$voteshare_last[i] <- match$voteshare
      data6$run_and_win_last[i] <- ifelse(match$voteshare>=0.5,1,0)
    }
    if(!(data6$cand_name[i] %in% prevelec$cand_name)){# if cand didn't run in last elec
      data6$run_last[i] <- 0
      data6$voteshare_last[i] <- NA
      data6$run_and_win_last[i] <- 0
    }
  }
}
summary(data6$voteshare_last)
summary(data6$run_last)
summary(data6$run_and_win_last)

(cct.placebo.last <- with(data6, rdrobust(y=voteshare_last, x=voteshare, c=0.5,
                                          bwselect="CCT", all=TRUE, kernel='uni')))


(cct.placebo.last.r <- with(data6, rdrobust(y=run_last, x=voteshare, c=0.5,
                                            bwselect="CCT", all=TRUE, kernel='uni')))

(cct.placebo.last.rw <- with(data6, rdrobust(y=run_and_win_last, x=voteshare, c=0.5,
                                             bwselect="CCT", all=TRUE, kernel='uni')))

####################
##### Table 1: #####
####################
#### Placebo Tests table:
rd.export5(list(cct.placebo.last, cct.placebo.last.rw, cct.placebo.last.r))
####


##############################
#### Descriptive analysis ####
##############################


data$mayor_incumbent <- NA
data$runnerup_incumbent <- NA
prev_elec <- NULL
for(i in 1:nrow(data)){
  prev_elec <- data[match(data$elec_index[i],data$elec_index_next),]
  if(nrow(prev_elec)>0){
    if(!is.na(data$mayor_name_final[i])){
      if(data$mayor_name_final[i] %in% prev_elec[,c("mayor_name_final")]){
        data$mayor_incumbent[i] <- 1 # marker for incumbents
      }
      if(!(data$mayor_name_final[i] %in% prev_elec[,c("mayor_name_final")])){
        data$mayor_incumbent[i] <- 0
      }
    }
    if(!is.na(data$runnerup_name_final[i])){
      if(data$runnerup_name_final[i] %in% prev_elec[,c("mayor_name_final")]){
        data$runnerup_incumbent[i] <- 1 # marker for incumbents
      }
      if(!(data$runnerup_name_final[i] %in% prev_elec[,c("mayor_name_final")])){
        data$runnerup_incumbent[i] <- 0
      }
    }
  }
  if(nrow(prev_elec)<1){
    data$mayor_incumbent[i] <- NA
    data$runnerup_incumbent[i] <- NA
  }
}

#####################
##### Figure 3: #####
#####################
pdf("voteshare_desc2.pdf")
par(mar=c(4,4,0,0))
plot(jitter((data$mayor_incumbent[!is.na(data$mayor_voteshare)]),2),data$mayor_voteshare[!is.na(data$mayor_voteshare)], col=ifelse(data$mayor_incumbent[!is.na(data$mayor_voteshare)]==1,"blue","black"), pch=ifelse(data$mayor_incumbent[!is.na(data$mayor_voteshare)]==1,1,17),axes=F,xlab="",ylab="Voteshare",main="",ylim=c(0,1))
points(jitter((data$runnerup_incumbent[!is.na(data$runnerup_voteshare)]),2),data$runnerup_voteshare[!is.na(data$runnerup_voteshare)], col=ifelse(data$runnerup_incumbent[!is.na(data$runnerup_voteshare)]==1,"blue","black"), pch=ifelse(data$runnerup_incumbent[!is.na(data$runnerup_voteshare)]==1,1,17))
axis(1,at=c(0,1),labels=c("Non-Incumbents","Incumbents"),tick = F,cex.axis=1.5)
axis(2,at=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),labels=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),cex.axis=1.5)
dev.off()
#####################

# Average percent of vote:
mean(c(data$mayor_voteshare[data$mayor_incumbent==1],data$runnerup_voteshare[data$runnerup_incumbent==1]),na.rm=T) # 65.42%
mean(c(data$mayor_voteshare[data$mayor_incumbent==0],data$runnerup_voteshare[data$runnerup_incumbent==0]),na.rm=T) # 44.10%
# Average vote margin:
mean((data$mayor_votes_final-data$runnerup_votes_final)/(data$mayor_votes_final+data$runnerup_votes_final),na.rm=T) # 31.81%
# Difference in average probability of running and winning, descriptively:
mean(data6$run_and_win_next[data6$voteshare>=0.5],na.rm=T) # 50.08
mean(data6$run_and_win_next[data6$voteshare<0.5],na.rm=T) # 3.40%

mean(data6$run_and_win_next[which(data6$run_next==1 & data6$voteshare>=0.5)],na.rm=T) # 80.83% chance of winning if they run again
mean(data6$run_and_win_next[which(data6$run_next==1 & data6$voteshare<0.5)],na.rm=T) # 21.67% chance of winning if they run again


########################
##### RDD analysis #####
########################
set.seed <- 02139

(cct.run <- with(data6, rdrobust(y=run_next , x=voteshare , c=0.5,
                                      bwselect="CCT", all=TRUE, kernel='uni')))


(cct.run.win <- with(data6, rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                bwselect="CCT", all=TRUE, kernel='uni')))


####################
##### Table 2: #####
####################
rd.export5(list(cct.run,cct.run.win))
####################


####################
##### Figure 4: ####
####################
## Running:
width <- .005
data6 <- mutate(data6, bin=cut(voteshare, breaks=seq(0, 1, width)))
bin2.df <- data.frame(bin.mean=tapply(data6$run_next[which(data6$voteshare>=(0.5-cct.run$h) & data6$voteshare<=(0.5+cct.run$h))],
                                      data6$bin[which(data6$voteshare>=(0.5-cct.run$h) & data6$voteshare<=(0.5+cct.run$h))], mean, na.rm=TRUE),
                      mid=seq(0 + width/2, 1 - width/2, width),
                      n=tapply(data6$run_next[which(data6$voteshare>=(0.5-cct.run$h) & data6$voteshare<=(0.5+cct.run$h))], data6$bin[which(data6$voteshare>=(0.5-cct.run$h) & data6$voteshare<=(0.5+cct.run$h))],
                               length))
tri <- function (x, h, c=0.5) pmax(0, 1 - abs((x - c) / h))

pdf("mayoral_rd_run2.pdf", width=8, height=6)
(ggplot(data6)
  + geom_point(aes(x = voteshare, y = run_next), shape="|", alpha=.5)
  + geom_smooth(data=subset(data6, voteshare > 0.5),
                aes(x = voteshare, y = run_next,
                    weight=tri(voteshare, cct.run$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5)
  + geom_smooth(data=subset(data6, voteshare < 0.5),
                aes(x = voteshare, y = run_next,
                    weight=tri(voteshare, cct.run$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5)
  + geom_point(data=bin2.df, aes(x=mid, y=bin.mean, size=n), pch=1)
  + geom_vline(xintercept=0)
  + geom_hline(yintercept=0, lty='dotted')
  + coord_cartesian(ylim = c(-0.05,1.05))
  + scale_x_continuous(breaks=seq(0,1,0.02),
                       limits=c(0.5-cct.run$h,
                                0.5+cct.run$h))
  
  + labs(size = '# Obs. in\n2-unit Bin')
  + theme_bw()
  + labs(x = "Vote share, time t",
         y = "Prob. Running, time t+1")
  #  + ggtitle("All Elections, Running")
  + annotate('text', x = 0.445, y = 0.75, 0.05, hjust=0, parse=TRUE,size=6,
             label = paste('hat(tau)==',
                           round(cct.run$coef['Conventional', ], 2),
                           '~(list(', round(cct.run$ci['Robust', 1], 2),
                           ',', round(cct.run$ci['Robust', 2], 2),
                           '))'))
  + theme(axis.text = element_text(size=15),axis.title = element_text(size=15))
  
)
dev.off()

## Running + Winning:
width <- .005
data6 <- mutate(data6, bin=cut(voteshare, breaks=seq(0, 1, width)))
bin2.df <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win$h) & data6$voteshare<=(0.5+cct.run.win$h))],
                                      data6$bin[which(data6$voteshare>=(0.5-cct.run.win$h) & data6$voteshare<=(0.5+cct.run.win$h))], mean, na.rm=TRUE),
                      mid=seq(0 + width/2, 1 - width/2, width),
                      n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win$h) & data6$voteshare<=(0.5+cct.run.win$h))], data6$bin[which(data6$voteshare>=(0.5-cct.run.win$h) & data6$voteshare<=(0.5+cct.run.win$h))],
                               length))
tri <- function (x, h, c=0.5) pmax(0, 1 - abs((x - c) / h))

pdf("mayoral_rd_runwin2.pdf", width=8, height=6)
(ggplot(data6)
  + geom_point(aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5)
  + geom_smooth(data=subset(data6, voteshare > 0.5),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5)
  + geom_smooth(data=subset(data6, voteshare < 0.5),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5)
  + geom_point(data=bin2.df, aes(x=mid, y=bin.mean, size=n), pch=1)
  + geom_vline(xintercept=0)
  + geom_hline(yintercept=0, lty='dotted')
  + coord_cartesian(ylim = c(-0.05,1.05))
  + scale_x_continuous(breaks=seq(0,1,0.02),
                       limits=c(0.5-cct.run.win$h,
                                0.5+cct.run.win$h))
  
  + labs(size = '# Obs. in\n2-unit Bin')
  + theme_bw()
  + labs(x = "Vote share, time t",
         y = "Prob. Running + Winning, time t+1")
  #  + ggtitle("All Elections, Running + Winning")
  + annotate('text', x = 0.445, y = 0.6, 0.05, hjust=0, parse=TRUE,size=6,
             label = paste('hat(tau)==',
                           round(cct.run.win$coef['Conventional', ], 2),
                           '~(list(', round(cct.run.win$ci['Robust', 1], 2),
                           ',', round(cct.run.win$ci['Robust', 2], 2),
                           '))'))
  + theme(axis.text = element_text(size=15),axis.title = element_text(size=15))
  
)
dev.off()
####################



#########################################
##### The Impact of Election Timing #####
#########################################

## Turnout differences across election timing:

data$turnout[data$turnout==0] <- NA
mean(data$turnout[data$concurrent==1],na.rm=T) # 37346.25
mean(data$turnout[data$concurrent==0],na.rm=T) # 30043.03, but confounded by population size

data$turnout_perc <- data$turnout/data$population_est
mean(data$turnout_perc[data$concurrent==1],na.rm=T) # 28.15%
mean(data$turnout_perc[data$concurrent==0],na.rm=T) # 17.90%


######
## Main results by Election Timing:
######

cct.run.on <- with(data6[data6$concurrent_next==1,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                             bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.off <- with(data6[data6$concurrent_next==0,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                              bwselect="CCT", all=TRUE, kernel='uni'))

(Z_run_timing = (cct.run.on$coef[1,1] - cct.run.off$coef[1,1])/sqrt(cct.run.on$se[3,1]^2 + cct.run.off$se[3,1]^2))
SE_run_timing <- sqrt(cct.run.on$se[1,1]^2 + cct.run.off$se[1,1]^2)
(Pvalue_run_timing <- 2*pnorm(-abs(Z_run_timing))) 


cct.run.win.on <- with(data6[data6$concurrent_next==1,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                    bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.off <- with(data6[data6$concurrent_next==0,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                             bwselect="CCT", all=TRUE, kernel='uni'))


####################
##### Table 3: #####
####################
rd.export5(list(cct.run.on, cct.run.off, cct.run.win.on, cct.run.win.off))
####################
####################

(Z_timing = (cct.run.win.on$coef[1,1] - cct.run.win.off$coef[1,1])/sqrt(cct.run.win.on$se[3,1]^2 + cct.run.win.off$se[3,1]^2))
SE_timing <- sqrt(cct.run.win.on$se[1,1]^2 + cct.run.win.off$se[1,1]^2)
(Pvalue_timing <- 2*pnorm(-abs(Z_timing))) # 0.0098

# check with fully interacted model:
summary(lm(run_and_win_next ~ voteshare_adj + (voteshare_adj>0) + voteshare_adj*(voteshare_adj>0) 
           + concurrent_next + voteshare_adj*concurrent_next + (voteshare_adj>0)*concurrent_next + voteshare_adj*(voteshare_adj>0)*concurrent_next,
           data=data6[which(data6$voteshare_adj<cct.run.win.off$h & data6$voteshare_adj>-cct.run.win.off$h),]))
# interaction significant

summary(lm(run_and_win_next ~ voteshare_adj + (voteshare_adj>0) + voteshare_adj*(voteshare_adj>0) 
           + concurrent_next + voteshare_adj*concurrent_next + (voteshare_adj>0)*concurrent_next + voteshare_adj*(voteshare_adj>0)*concurrent_next,
           data=data6[which(data6$voteshare_adj<cct.run.win.on$h & data6$voteshare_adj>-cct.run.win.on$h),]))
# interaction significant



####################
##### Figure 5: ####
####################

width <- .005
data6 <- mutate(data6, bin=cut(voteshare, breaks=seq(0, 1, width)))
bin2.df <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.on$h) & data6$voteshare<=(0.5+cct.run.win.on$h) & data6$concurrent_next==1)],
                                      data6$bin[which(data6$voteshare>=(0.5-cct.run.win.on$h) & data6$voteshare<=(0.5+cct.run.win.on$h) & data6$concurrent_next==1)], mean, na.rm=TRUE),
                      mid=seq(0 + width/2, 1 - width/2, width),
                      n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.on$h) & data6$voteshare<=(0.5+cct.run.win.on$h) & data6$concurrent_next==1)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.on$h) & data6$voteshare<=(0.5+cct.run.win.on$h) & data6$concurrent_next==1)],
                               length))
bin2.dfb <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.on$h) & data6$voteshare<=(0.5+cct.run.win.on$h) & data6$concurrent_next==0)],
                                       data6$bin[which(data6$voteshare>=(0.5-cct.run.win.on$h) & data6$voteshare<=(0.5+cct.run.win.on$h) & data6$concurrent_next==0)], mean, na.rm=TRUE),
                       mid=seq(0 + width/2, 1 - width/2, width),
                       n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.on$h) & data6$voteshare<=(0.5+cct.run.win.on$h) & data6$concurrent_next==0)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.on$h) & data6$voteshare<=(0.5+cct.run.win.on$h) & data6$concurrent_next==0)],
                                length))
tri <- function (x, h, c=0.5) pmax(0, 1 - abs((x - c) / h))

pdf("mayoral_rd_runwin_timing.pdf", width=8, height=6)
(ggplot(data6)
  + geom_point(data=subset(data6,concurrent_next == 0),
               aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5)
  + geom_point(data=subset(data6,concurrent_next == 1),
               aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5,col="red")
  + geom_smooth(data=subset(data6, voteshare > 0.5 & concurrent_next == 0),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.off$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
  + geom_smooth(data=subset(data6, voteshare > 0.5 & concurrent_next == 1),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.on$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
  + geom_smooth(data=subset(data6, voteshare < 0.5 & concurrent_next == 0),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.off$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
  + geom_smooth(data=subset(data6, voteshare < 0.5 & concurrent_next == 1),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.on$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
  + geom_point(data=bin2.dfb, aes(x=mid, y=bin.mean, size=n), pch=1)
  + geom_point(data=bin2.df, aes(x=mid, y=bin.mean, size=n), pch=1,col="red")
  + geom_vline(xintercept=0)
  + geom_hline(yintercept=0, lty='dotted')
  + coord_cartesian(ylim = c(-0.05,1.05))
  + scale_x_continuous(breaks=seq(0,1,0.01),
                       limits=c(0.5-cct.run.win.on$h,
                                0.5+cct.run.win.on$h))
  
  + labs(size = '# Obs. in\n2-unit Bin', color="Election Timing")
  + theme_bw()
  + labs(x = "Vote share, time t",
         y = "Prob. Running + Winning, time t+1")
  #  + ggtitle("Running + Winning by Election Timing")
  + annotate('text', x = 0.45, y = 0.33, 0.05, hjust=0, parse=TRUE, col="black",size=6,
             label = paste('hat(tau)==',
                           round(cct.run.win.off$coef['Conventional', ], 2),
                           '~(list(', round(cct.run.win.off$ci['Robust', 1], 2),
                           ',', round(cct.run.win.off$ci['Robust', 2], 2),
                           '))'))
  + annotate('text', x = 0.45, y = 0.74, 0.05, hjust=0, parse=TRUE, col="red",size=6,
             label = paste('hat(tau)==',
                           round(cct.run.win.on$coef['Conventional', ], 2),
                           '~(list(', round(cct.run.win.on$ci['Robust', 1], 2),
                           ',', round(cct.run.win.on$ci['Robust', 2], 2),
                           '))'))
  + annotate('text', x = 0.45, y = 0.382, 0.05, hjust=0, parse=TRUE,
             label = "Off-cycle",col="black",size=6)
  + annotate('text', x = 0.45, y = 0.792, 0.05, hjust=0, parse=TRUE,
             label = "On-cycle",col="red",size=6)
  + theme(axis.text = element_text(size=15),axis.title = element_text(size=15))
  
)
dev.off()
####################


#################
### Changes in Election Timing:
#################
# create subgroups:
allelecs <- data6[!duplicated(data6$elec_index),]
switchers <- ddply(allelecs, ~ fips, summarize,
                   var_timing = var(concurrent,na.rm=T))
switchers <- switchers[switchers$var_timing>0,] # leaves 216 cities
switchers <- switchers[!is.na(switchers$fips),] # leaves 76 cities
switchcities <- allelecs[which(allelecs$fips %in% switchers$fips),] # 1075 elections in those cities



##############################
##### Inc. Adv. by group #####
##############################

data6$switcher <- ifelse(data6$fips %in% switchers$fips,1,0) # 1092 elections in those cities
table(data6$switcher) # 1969 cand-elecs that were switchers

cct.switchers.on <-  with(data6[which(data6$concurrent_next==1 & data6$switcher==1),],
                          rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                   bwselect="CCT", all=TRUE, kernel='uni'))

cct.switchers.off <-  with(data6[which(data6$concurrent_next==0 & data6$switcher==1),],
                           rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                    bwselect="CCT", all=TRUE, kernel='uni'))

cct.constant.on <-  with(data6[which(data6$concurrent_next==1 & data6$switcher==0),],
                         rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                  bwselect="CCT", all=TRUE, kernel='uni'))

cct.constant.off <-  with(data6[which(data6$concurrent_next==0 & data6$switcher==0),],
                          rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                   bwselect="CCT", all=TRUE, kernel='uni'))


cct.constant.on$coef[1,1] - cct.constant.off$coef[1,1] # 45.72% difference
cct.switchers.on$coef[1,1] - cct.switchers.off$coef[1,1] # 26.55% difference

(Z_timing.switchers = (cct.switchers.on$coef[1,1] - cct.switchers.off$coef[1,1])/sqrt(cct.switchers.on$se[3,1]^2 + cct.switchers.off$se[3,1]^2))
SE_timing.switchers <- sqrt(cct.switchers.on$se[1,1]^2 + cct.switchers.off$se[1,1]^2)
(Pvalue_timing.switchers <- 2*pnorm(-abs(Z_timing.switchers))) # 0.052

(Z_timing.constant = (cct.constant.on$coef[1,1] - cct.constant.off$coef[1,1])/sqrt(cct.constant.on$se[3,1]^2 + cct.constant.off$se[3,1]^2))
SE_timing.constant <- sqrt(cct.constant.on$se[1,1]^2 + cct.constant.off$se[1,1]^2)
(Pvalue_timing.constant <- 2*pnorm(-abs(Z_timing.constant))) # 0.000

### make into table:
tab <- cbind(rbind("Non-switchers","","","Switchers","",""),
             rbind(round(cct.constant.on$coef[1,1],2),paste0("(",round(cct.constant.on$ci[3,1],2),", ",round(cct.constant.on$ci[3,2],2),")"),paste0("n = ",cct.constant.on$N),
                   round(cct.switchers.on$coef[1,1],2),paste0("(",round(cct.switchers.on$ci[3,1],2),", ",round(cct.switchers.on$ci[3,2],2),")"),paste0("n = ",cct.switchers.on$N)),
             rbind(round(cct.constant.off$coef[1,1],2),paste0("(",round(cct.constant.off$ci[3,1],2),", ",round(cct.constant.off$ci[3,2],2),")"),paste0("n = ",cct.constant.off$N),
                   round(cct.switchers.off$coef[1,1],2),paste0("(",round(cct.switchers.off$ci[3,1],2),", ",round(cct.switchers.off$ci[3,2],2),")"),paste0("n = ",cct.switchers.off$N)),
             rbind(round(cct.constant.on$coef[1,1]-cct.constant.off$coef[1,1],2),"","",round(cct.switchers.on$coef[1,1]-cct.switchers.off$coef[1,1],2),"",""),
             rbind(round(Pvalue_timing.constant,2),"","",round(Pvalue_timing.switchers,2),"","")
)
colnames(tab) <- c("Subset","On-cycle Inc. Adv.","Off-cycle Inc. Adv.","Difference","p-value of difference")

####################
##### Table 4: #####
print(xtable(tab),include.rownames = F,booktabs=T)
####################

inc_df_switchers <- data.frame(subset=c("cities that stay on-cycle",
                                        "cities that stay off-cycle",
                                        "cities that switch timing",
                                        "cities that switch timing"),
                               timing=c("on","off","on","off"),
                               inc_adv=c(cct.constant.on$coef[1,1],
                                         cct.constant.off$coef[1,1],
                                         cct.switchers.on$coef[1,1],
                                         cct.switchers.off$coef[1,1]),
                               se=c(cct.constant.on$se[3,1],
                                    cct.constant.off$se[3,1],
                                    cct.switchers.on$se[3,1],
                                    cct.switchers.off$se[3,1]),
                               ci.lo=c(cct.constant.on$ci[3,1],
                                       cct.constant.off$ci[3,1],
                                       cct.switchers.on$ci[3,1],
                                       cct.switchers.off$ci[3,1]),
                               ci.hi=c(cct.constant.on$ci[3,2],
                                       cct.constant.off$ci[3,2],
                                       cct.switchers.on$ci[3,2],
                                       cct.switchers.off$ci[3,2]),
                               n=c(cct.constant.on$N,
                                   cct.constant.off$N,
                                   cct.switchers.on$N,
                                   cct.switchers.off$N))

####################
##### Figure 6: ####
####################
pdf("diffs_timing_grouped.pdf",height=5,width=7)
par(mar=c(4,5,1,1))
plot("",xlim=c(0.75,2.25),ylim=c(0,1.2),axes=F,
     xlab="Timing",
     ylab="Incumbency Advantage\n(Probability Running + Winning)")
segments(x0=0.5,x1=2.5,y0=seq(0.0,1.2,0.1),y1=seq(0.0,1.2,0.1),col="lightgrey")
points(c(1,2,1.2,1.8),inc_df_switchers[,3],pch=c(16,16,17,17),col=c("red","black","blue","blue"))
segments(x0=c(1,2,1.2,1.8),x1=c(1,2,1.2,1.8),y0 = inc_df_switchers$ci.lo, y1 = inc_df_switchers$ci.hi,col=c("red","black","blue","blue"))
segments(x0=1.2,x1=1.8,y0 = inc_df_switchers[3,3], y1 = inc_df_switchers[4,3],col="blue")
axis(2,at=seq(0,1.0,0.2),labels = seq(0,1.0,0.2),tick = T)
axis(1,at=c(1.1,1.9),labels = c("On-cycle","Off-cycle"),tick = T)
text(1.5,0.2,"Cities that\nswitch timing",col="blue")
text(1.05,0.86,"Cities that\ndo not switch\ntiming",col="red",adj=0)
text(2.05,0.365,"Cities that\ndo not switch\ntiming",col="black",adj=0,xpd=T)
dev.off()
####################


data$switcher <- ifelse(data$fips %in% switchers$fips,1,0) 
table(data$switcher) # 1092 elections in ...
length(unique(data$fips[which(data$switcher==1)])) # ... 76 cities

turnout_summ2 <- ddply(data[!is.na(data$concurrent),],~ concurrent + switcher,summarize,
                       meanturnout = mean(turnout_perc,na.rm=T),
                       SDturnout = sd(turnout_perc,na.rm=T),
                       n = length(turnout_perc))
turnout_summ2$ci.lo <- turnout_summ2$meanturnout-turnout_summ2$SDturnout
turnout_summ2$ci.hi <- turnout_summ2$meanturnout+turnout_summ2$SDturnout
t.test(data$turnout_perc[data$concurrent==1 & data$switcher==0],
       data$turnout_perc[data$concurrent==0 & data$switcher==0]) # significant difference
t.test(data$turnout_perc[data$concurrent==1 & data$switcher==1],
       data$turnout_perc[data$concurrent==0 & data$switcher==1]) # significant difference



data$margin <- (data$mayor_votes_final-data$runnerup_votes_final)
data$margin_perc <- (data$mayor_votes_final-data$runnerup_votes_final)/(data$mayor_votes_final+data$runnerup_votes_final)

mean(data$margin[data$concurrent==1]/data$population_est[data$concurrent==1],na.rm=T) # 9.1%
mean(data$margin[data$concurrent==0]/data$population_est[data$concurrent==0],na.rm=T) # 4.68% margin -  4.5pp larger in on-cycle
t.test(data$margin[data$concurrent==1]/data$population_est[data$concurrent==1],data$margin[data$concurrent==0]/data$population_est[data$concurrent==0]) # sig different (p<0.0001)

mean(data$margin[data$concurrent==1],na.rm=T)
mean(data$margin[data$concurrent==0],na.rm=T) # large different in number of votes
t.test(data$margin[data$concurrent==1],data$margin[data$concurrent==0]) # sig different (p<0.000)

########
## Mechanisms:
########
# Experiment data analysis in separate file




#####################
# # # # # # # # # # #
# # # Appendix # # #
# # # # # # # # # # #
#####################


##### Appendix A
# Robustness checks:

## need to sample one candidate from every election:
set.seed <- 02139
sample_data <- NULL
for(i in 1:length(unique(data6$elec_index))){
  thiselec <- data6[which(data6$elec_index==data6$elec_index[i]),]
  thiscand <- thiselec[sample(1:nrow(thiselec),1),]
  sample_data <- rbind(sample_data,thiscand)
}
dim(sample_data) # 8974 candidate-elections


# formal McCrary test:
width <- .005
sample_data <- mutate(sample_data, bin=cut(voteshare_adj, breaks=seq(-0.5, 0.5, width)))
bin2.df <- data.frame(bin.mean=tapply(sample_data$voteshare_adj,
                                      sample_data$bin, mean, na.rm=TRUE),
                      mid=seq(-0.5 + width/2, 0.5 - width/2, width),
                      n=tapply(sample_data$voteshare_adj, sample_data$bin,
                               length))
mc.rd <- lm(n ~ mid * (mid>=0), data=bin2.df[which(bin2.df$mid>=(-cct.run.win$h) & bin2.df$mid<=(cct.run.win$h)),])
summary(mc.rd) 


### Density histogram with sampled data:

mccrary_bw <- ggplot(sample_data[which(sample_data$voteshare_adj>-0.5 & sample_data$voteshare_adj<0.5),],aes(x=voteshare_adj)) +
  geom_histogram(aes(y=..density..), alpha=0.8, binwidth=0.03, col="white") +
  geom_density(col="cornflowerblue",lwd=2) + 
  labs(x="Vote Margin",y="Density") +
  xlim(c(-0.55,0.55)) + 
  ylim(c(0,3)) +
  scale_x_continuous(breaks=seq(-0.5,0.5,0.1),
                     limits=c(-0.55,0.55)) +
  theme(text = element_text(size=20)) + 
  theme_bw()

#######################
##### Figure A-1: #####
#######################
ggsave(mccrary_bw,filename = "mccrary_sample_bw.pdf",height=7,width=14)
#######################




### Appendix B: Bandwidth checks

data6$bandwidth <- abs(data6$voteshare-0.5)

test_bandwidths <- matrix(nrow=50, ncol=3)
colnames(test_bandwidths) <- c("bandwidth", "coef", "se")
test_bandwidths <- as.data.frame(test_bandwidths)
test_bandwidths$bandwidth <- c(seq(.01,.09,.01),cct.run.win$h,seq(.11,.5,.01))

for(i in 1:length(test_bandwidths$bandwidth)){
  reg <- with(data6[which(data6$bandwidth < test_bandwidths$bandwidth[i]),], lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>0)))
  test_bandwidths$coef[i] <- coef(reg)[3]
  test_bandwidths$se[i] <- coef(summary(reg))[3,2]
}

test_bandwidths$ci <- paste0("(",round(test_bandwidths$coef + qnorm(0.025)*test_bandwidths$se,2),", ",round(test_bandwidths$coef + qnorm(0.975)*test_bandwidths$se,2),")")

######################
##### Table B-1: #####
######################
print(xtable(test_bandwidths[,c("bandwidth","coef","ci")]),include.rownames=F)
######################


test_bandwidths$color <- "red"
test_bandwidths$color[test_bandwidths$bandwidth==cct.run.win$h] <- "blue"

title <- paste0("Incumbency Advantage, Probability of Running + Winning,\nSensitivity to Different Bandwidths\n")

#######################
##### Figure B-2: #####
#######################
pdf("test_bandwidths_all.pdf",width=10,height=7) # for paper
ggplot(data = test_bandwidths, aes(y = coef, x = bandwidth*100)) +
  geom_point( colour = test_bandwidths$color, size = 3)+
  geom_errorbar(aes(x = bandwidth*100, 
                    ymin = coef - 1.96*se, ymax = coef + 1.96*se, height=.1),
                colour = test_bandwidths$color) +
  xlab("Bandwidth (optimal bandwidth to minimize bias plotted in blue)") + ylab("Incumbency Advantage\n(Probability of Running + Winning)") +
  expand_limits(y=c(0,1.0)) +
  geom_vline(xintercept = 0,size=.5,colour="blue",linetype="dotted") +
  theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black")) +  
  theme_bw() +
  theme(plot.title = element_text(lineheight=.8, face="bold"),axis.text = element_text(size=15),axis.title = element_text(size=15))
dev.off()
#######################





## Incumbency advantage by time period:

# before 1972
cct.early <- with(data6[data6$YearData<1972,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                        bwselect="CCT", all=TRUE, kernel='uni'))
## 1972-1991
cct.mid <- with(data6[data6$YearData>1971 & data6$YearData<1992,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                            bwselect="CCT", all=TRUE, kernel='uni'))
## 1992-2014
cct.late <- with(data6[data6$YearData>1991,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                       bwselect="CCT", all=TRUE, kernel='uni'))


######################
##### Table B-2: #####
######################
rd.export5(list(cct.early,cct.mid,cct.late))
######################

# significance tests:
(cct.late$coef[1,1] - cct.early$coef[1,1]) # 0.028
(Z_time1 = (cct.late$coef[1,1] - cct.early$coef[1,1])/sqrt(cct.late$se[3,1]^2 + cct.early$se[3,1]^2))
(Pvalue_time1 <- 2*pnorm(-abs(Z_time1))) # 0.6299
(cct.mid$coef[1,1] - cct.early$coef[1,1]) # 0.076
(Z_time2 = (cct.mid$coef[1,1] - cct.early$coef[1,1])/sqrt(cct.mid$se[3,1]^2 + cct.early$se[3,1]^2))
(Pvalue_time2 <- 2*pnorm(-abs(Z_time2))) # 0.177
(cct.late$coef[1,1] - cct.mid$coef[1,1])
(Z_time3 = (cct.late$coef[1,1] - cct.mid$coef[1,1])/sqrt(cct.late$se[3,1]^2 + cct.mid$se[3,1]^2))
(Pvalue_time3 <- 2*pnorm(-abs(Z_time3))) # 0.362


## by city size:

# small cities
cct.small <- with(data6[which(data6$population_est<75000),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                      bwselect="CCT", all=TRUE, kernel='uni'))

## medium cities
cct.medium <- with(data6[which(data6$population_est<170000 & data6$population_est>=75000),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                                                      bwselect="CCT", all=TRUE, kernel='uni'))

## large cities (>170k population)
cct.large <- with(data6[data6$population_est>=170000,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                 bwselect="CCT", all=TRUE, kernel='uni'))

######################
##### Table B-3: #####
######################
rd.export5(list(cct.small,cct.medium,cct.large))
######################





# significance tests:
(cct.large$coef[1,1] - cct.small$coef[1,1]) # 0.055
Z_size1 <- (cct.small$coef[1,1] - cct.large$coef[1,1])/sqrt(cct.small$se[3,1]^2 + cct.large$se[3,1]^2)
SE_size1 <- sqrt(cct.small$se[1,1]^2 + cct.large$se[1,1]^2)
(Pvalue_size1 <- 2*pnorm(-abs(Z_size1))) # 0.387
(cct.medium$coef[1,1] - cct.small$coef[1,1]) # 0.013
Z_size3 <- (cct.medium$coef[1,1] - cct.small$coef[1,1])/sqrt(cct.medium$se[3,1]^2 + cct.small$se[3,1]^2)
SE_size3 <- sqrt(cct.medium$se[1,1]^2 + cct.small$se[1,1]^2)
(Pvalue_size3 <- 2*pnorm(-abs(Z_size3))) # 0.843
(cct.large$coef[1,1] - cct.medium$coef[1,1]) # 0.042
Z_size2 <- (cct.medium$coef[1,1] - cct.large$coef[1,1])/sqrt(cct.medium$se[3,1]^2 + cct.large$se[3,1]^2)
SE_size2 <- sqrt(cct.medium$se[1,1]^2 + cct.large$se[1,1]^2)
(Pvalue_size2 <- 2*pnorm(-abs(Z_size2))) # 0.559



### Run tests of bandwidth size by election timing:
data6$bandwidth <- abs(data6$voteshare-0.5)

test_bandwidths <- matrix(nrow=52, ncol=6)
colnames(test_bandwidths) <- c("bandwidth", "coef_on", "se_on","coef_off","se_off","p_value")
test_bandwidths <- as.data.frame(test_bandwidths)
test_bandwidths$bandwidth <- c(seq(.01,.05,.01),cct.run.win.on$h,seq(.06,.1,.01),cct.run.win.off$h,seq(0.11,0.5,.01))

for(i in 1:length(test_bandwidths$bandwidth)){
  reg.on <- with(data6[which(data6$bandwidth < test_bandwidths$bandwidth[i] & data6$concurrent_next==1),], lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>0)))
  test_bandwidths$coef_on[i] <- coef(reg.on)[3]
  test_bandwidths$se_on[i] <- coef(summary(reg.on))[3,2]
  reg.off <- with(data6[which(data6$bandwidth < test_bandwidths$bandwidth[i] & data6$concurrent_next==0),], lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>0)))
  test_bandwidths$coef_off[i] <- coef(reg.off)[3]
  test_bandwidths$se_off[i] <- coef(summary(reg.off))[3,2]
  Z_timing <- (reg.on$coef[3] - reg.off$coef[3])/sqrt(summary(reg.on)$coefficients[3,2]^2 + summary(reg.off)$coefficients[3,2]^2)
  test_bandwidths$p_value[i] <- 2*pnorm(-abs(Z_timing))
}

test_bandwidths$color_on <- "black"
test_bandwidths$color_on[test_bandwidths$bandwidth==cct.run.win.on$h] <- "blue"
test_bandwidths$color_off <- "red"
test_bandwidths$color_off[test_bandwidths$bandwidth==cct.run.win.off$h] <- "blue"
test_bandwidths$coef_on[test_bandwidths$bandwidth==0.06] <- test_bandwidths$se_on[test_bandwidths$bandwidth==0.06] <- NA
test_bandwidths$coef_off[test_bandwidths$bandwidth==0.11] <- test_bandwidths$se_off[test_bandwidths$bandwidth==0.11] <- NA
test_bandwidths$coef_on[test_bandwidths$bandwidth==cct.run.win.off$h] <- test_bandwidths$se_on[test_bandwidths$bandwidth==cct.run.win.off$h] <- NA
test_bandwidths$coef_off[test_bandwidths$bandwidth==cct.run.win.on$h] <- test_bandwidths$se_off[test_bandwidths$bandwidth==cct.run.win.on$h] <- NA


#######################
##### Figure B-3: #####
#######################
pdf(paste0("test_bandwidths_timing2",".pdf"),width=10,height=7) # for paper
ggplot(data = test_bandwidths, aes(x = bandwidth*100)) +
  geom_ribbon(data=subset(test_bandwidths,!is.na(se_on)),aes(x = bandwidth*100, 
                                                             ymin = coef_on - 1.96*se_on, ymax = coef_on + 1.96*se_on, height=.1),
              colour = "gray",alpha=0.1) +
  geom_ribbon(data=subset(test_bandwidths,!is.na(se_off)), aes(x = bandwidth*100, 
                                                               ymin = coef_off - 1.96*se_off, ymax = coef_off + 1.96*se_off, height=.1),
              colour = "gray", alpha=0.1) +
  geom_point(aes(y=coef_on),colour = test_bandwidths$color_on, size = 3,shape=16) +
  geom_point(aes(y=coef_off),colour = test_bandwidths$color_off, size = 3,shape=17) +
  xlab("Bandwidth (optimal bandwidths to minimize bias plotted in blue)") + ylab("Incumbency Advantage\n(Probability of Running + Winning)") +
  expand_limits(y=c(0,1.0)) +
  geom_vline(xintercept = 0,size=.5,colour="blue",linetype="dotted") +
  theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black")) +  
  #   ggtitle(title) +
  theme_bw() +
  theme(plot.title = element_text(lineheight=.8, face="bold"),axis.text = element_text(size=15),axis.title = element_text(size=15)) +
  annotate('text', x = 10, y = 0.75, hjust=0, col="black",
           label = paste0("On-cycle")) +
  annotate('text', x = 25, y = 0.25, hjust=0, col="red",
           label = paste0("Off-cycle"))
dev.off()
#######################


## Subgroup results
cct.run.on.pres <- with(data6[data6$concurrent_pres_next==1,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                                        bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.on.pres <- with(data6[data6$concurrent_pres_next==1,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                            bwselect="CCT", all=TRUE, kernel='uni'))


cct.run.on.mid <- with(data6[which(data6$concurrent_pres_next==0 & data6$concurrent_next==1),], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                                                                         bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.on.mid <- with(data6[which(data6$concurrent_pres_next==0 & data6$concurrent_next==1),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                                                             bwselect="CCT", all=TRUE, kernel='uni'))


#####################
##### Table B-4 #####
#####################
rd.export5(list(cct.run.win.on.mid,cct.run.win.on.pres,cct.run.win.off))

# significance tests:
(cct.run.win.on.pres$coef[1,1] - cct.run.win.off$coef[1,1])
Z_timing = (cct.run.win.on.pres$coef[1,1] - cct.run.win.off$coef[1,1])/sqrt(cct.run.win.on.pres$se[3,1]^2 + cct.run.win.off$se[3,1]^2)
SE_timing <- sqrt(cct.run.win.on.pres$se[1,1]^2 + cct.run.win.off$se[1,1]^2)
(Pvalue_timing <- 2*pnorm(-abs(Z_timing))) # 0.41

(cct.run.win.on.mid$coef[1,1] - cct.run.win.off$coef[1,1])
Z_timing = (cct.run.win.on.mid$coef[1,1] - cct.run.win.off$coef[1,1])/sqrt(cct.run.win.on.mid$se[3,1]^2 + cct.run.win.off$se[3,1]^2)
SE_timing <- sqrt(cct.run.win.on.mid$se[1,1]^2 + cct.run.win.off$se[1,1]^2)
(Pvalue_timing <- 2*pnorm(-abs(Z_timing))) # 0.002

################


######################
##### Appendix C #####
######################
# Party-level incumbency advantage
summary(data$demshare_next) # using this variable, which has been 0-coded if Dem party doesn't run in next election
data$demwin_next <- ifelse(data$demshare_next>=0.5,1,0)

(cct.partyIA.vs <- with(data,rdrobust(y=demshare_next , x=demshare , c=0.5)))

(cct.partyIA.pr <- with(data,rdrobust(y=demwin_next , x=demshare , c=0.5)))

#####################
##### Table C-5 #####
#####################
rd.export5(list(cct.partyIA.vs,cct.partyIA.pr))
#####################


## By timing:

(cct.partyIA.vs.on <- with(data[data$concurrent_next==1,], rdrobust(y=demshare_next , x=demshare , c=0.5)))

(cct.partyIA.vs.off <- with(data[data$concurrent_next==0,], rdrobust(y=demshare_next , x=demshare , c=0.5)))


(cct.partyIA.pr.on <- with(data[data$concurrent_next==1,],rdrobust(y=demwin_next , x=demshare , c=0.5)))

(cct.partyIA.pr.off <- with(data[data$concurrent_next==0,],rdrobust(y=demwin_next , x=demshare , c=0.5)))

#####################
##### Table C-6 #####
#####################
rd.export5(list(cct.partyIA.vs.on,cct.partyIA.vs.off,cct.partyIA.pr.on,cct.partyIA.pr.off))
#####################




######################
##### Appendix D #####
######################
# Additional institutional results

### Incumbency advantage by partisan ballots:
cct.run.p <- with(data6[data6$partisan==1,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                         bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.np <- with(data6[data6$partisan==0,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                          bwselect="CCT", all=TRUE, kernel='uni'))


cct.run.win.p <- with(data6[data6$partisan==1,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                             bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.np <- with(data6[data6$partisan==0,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                              bwselect="CCT", all=TRUE, kernel='uni'))

#####################
##### Table D-7: ####
#####################
rd.export5(list(cct.run.p,cct.run.np,cct.run.win.p,cct.run.win.np))
######################

(cct.run.win.p$coef[1,1] - cct.run.win.np$coef[1,1])
Z_pid = (cct.run.win.p$coef[1,1] - cct.run.win.np$coef[1,1])/sqrt(cct.run.win.p$se[3,1]^2 + cct.run.win.np$se[3,1]^2)
SE_pid <- sqrt(cct.run.win.p$se[1,1]^2 + cct.run.win.np$se[1,1]^2)
(Pvalue_pid <- 2*pnorm(-abs(Z_pid))) # 


######################
##### Figure D-4 #####
######################
width <- .005
data6 <- mutate(data6, bin=cut(voteshare, breaks=seq(0, 1, width)))
bin2.df <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.p$h) & data6$voteshare<=(0.5+cct.run.win.p$h) & data6$partisan==0)],
                                      data6$bin[which(data6$voteshare>=(0.5-cct.run.win.p$h) & data6$voteshare<=(0.5+cct.run.win.p$h) & data6$partisan==0)], mean, na.rm=TRUE),
                      mid=seq(0 + width/2, 1 - width/2, width),
                      n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.p$h) & data6$voteshare<=(0.5+cct.run.win.p$h) & data6$partisan==0)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.p$h) & data6$voteshare<=(0.5+cct.run.win.p$h) & data6$partisan==0)],
                               length))
bin2.dfb <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.p$h) & data6$voteshare<=(0.5+cct.run.win.p$h) & data6$partisan==1)],
                                       data6$bin[which(data6$voteshare>=(0.5-cct.run.win.p$h) & data6$voteshare<=(0.5+cct.run.win.p$h) & data6$partisan==1)], mean, na.rm=TRUE),
                       mid=seq(0 + width/2, 1 - width/2, width),
                       n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.p$h) & data6$voteshare<=(0.5+cct.run.win.p$h) & data6$partisan==1)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.p$h) & data6$voteshare<=(0.5+cct.run.win.p$h) & data6$partisan==1)],
                                length))
tri <- function (x, h, c=0.5) pmax(0, 1 - abs((x - c) / h))

pdf("mayoral_rd_runwin_pid.pdf", width=8, height=6)
(ggplot(data6)
  + geom_point(data=subset(data6,partisan == 0),
               aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5)
  + geom_point(data=subset(data6,partisan == 1),
               aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5,col="red")
  + geom_smooth(data=subset(data6, voteshare > 0.5 & partisan == 0),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.np$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
  + geom_smooth(data=subset(data6, voteshare > 0.5 & partisan == 1),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.p$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
  + geom_smooth(data=subset(data6, voteshare < 0.5 & partisan == 0),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.np$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
  + geom_smooth(data=subset(data6, voteshare < 0.5 & partisan == 1),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.p$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
  + geom_point(data=bin2.df, aes(x=mid, y=bin.mean, size=n), pch=1)
  + geom_point(data=bin2.dfb, aes(x=mid, y=bin.mean, size=n), pch=1,col="red")
  + geom_vline(xintercept=0)
  + geom_hline(yintercept=0, lty='dotted')
  + coord_cartesian(ylim = c(-0.05,1.05))
  + scale_x_continuous(breaks=seq(0,1,0.01),
                       limits=c(0.5-cct.run.win.p$h,
                                0.5+cct.run.win.p$h))
  
  + labs(size = '# Obs. in\n2-unit Bin')
  + theme_bw()
  + labs(x = "Vote share, time t",
         y = "Prob. Running + Winning, time t+1")
  #  + ggtitle("Running + Winning, Partisan vs. Nonpartisan Ballots")
  + annotate('text', x = 0.455, y = 0.33, 0.05, hjust=0, parse=TRUE, col="black",
             label = paste('hat(tau)==',
                           round(cct.run.win.np$coef['Conventional', ], 2),
                           '~(list(', round(cct.run.win.np$ci['Robust', 1], 2),
                           ',', round(cct.run.win.np$ci['Robust', 2], 2),
                           '))'))
  + annotate('text', x = 0.455, y = 0.6, 0.05, hjust=0, parse=TRUE, col="red",
             label = paste('hat(tau)==',
                           round(cct.run.win.p$coef['Conventional', ], 2),
                           '~(list(', round(cct.run.win.p$ci['Robust', 1], 2),
                           ',', round(cct.run.win.p$ci['Robust', 2], 2),
                           '))'))
  + annotate('text', x = 0.455, y = 0.38, 0.05, hjust=0, parse=TRUE,
             label = "Nonpartisan",col="black")
  + annotate('text', x = 0.455, y = 0.65, 0.05, hjust=0, parse=TRUE,
             label = "Partisan",col="red")
  
)
dev.off()
######################




### By form of government:
cct.run.may <- with(data6[data6$fog==1,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                      bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.man <- with(data6[data6$fog==2,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                       bwselect="CCT", all=TRUE, kernel='uni'))


cct.run.win.may <- with(data6[data6$fog==1,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                          bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.man <- with(data6[data6$fog==2,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                           bwselect="CCT", all=TRUE, kernel='uni'))
#####################
##### Table D-8: ####
#####################
rd.export5(list(cct.run.may,cct.run.man,cct.run.win.may,cct.run.win.man))
######################
(cct.run.win.may$coef[1,1] - cct.run.win.man$coef[1,1])
Z_fog = (cct.run.win.may$coef[1,1] - cct.run.win.man$coef[1,1])/sqrt(cct.run.win.may$se[3,1]^2 + cct.run.win.man$se[3,1]^2)
SE_fog <- sqrt(cct.run.win.may$se[1,1]^2 + cct.run.win.man$se[1,1]^2)
(Pvalue_fog <- 2*pnorm(-abs(Z_fog))) # 


######################
##### Figure D-5 #####
######################
### Form of Government:
width <- .005
data6 <- mutate(data6, bin=cut(voteshare, breaks=seq(0, 1, width)))
bin2.dfb <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.may$h) & data6$voteshare<=(0.5+cct.run.win.may$h) & data6$fog2==0)],
                                       data6$bin[which(data6$voteshare>=(0.5-cct.run.win.may$h) & data6$voteshare<=(0.5+cct.run.win.may$h) & data6$fog2==0)], mean, na.rm=TRUE),
                       mid=seq(0 + width/2, 1 - width/2, width),
                       n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.may$h) & data6$voteshare<=(0.5+cct.run.win.may$h) & data6$fog2==0)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.may$h) & data6$voteshare<=(0.5+cct.run.win.may$h) & data6$fog2==0)],
                                length))

bin2.df <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.man$h) & data6$voteshare<=(0.5+cct.run.win.man$h) & data6$fog2==1)],
                                      data6$bin[which(data6$voteshare>=(0.5-cct.run.win.man$h) & data6$voteshare<=(0.5+cct.run.win.man$h) & data6$fog2==1)], mean, na.rm=TRUE),
                      mid=seq(0 + width/2, 1 - width/2, width),
                      n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.man$h) & data6$voteshare<=(0.5+cct.run.win.man$h) & data6$fog2==1)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.man$h) & data6$voteshare<=(0.5+cct.run.win.man$h) & data6$fog2==1)],
                               length))
tri <- function (x, h, c=0.5) pmax(0, 1 - abs((x - c) / h))

pdf("mayoral_rd_runwin_fog.pdf", width=8, height=6)
(ggplot(data6)
  + geom_point(data=subset(data6,fog2 == 0), # city manager
               aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5)
  + geom_point(data=subset(data6,fog2 == 1), # strong mayor
               aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5,col="red")
  + geom_smooth(data=subset(data6, voteshare > 0.5 & fog2 == 0),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.man$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
  + geom_smooth(data=subset(data6, voteshare > 0.5 & fog2 == 1),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.may$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
  + geom_smooth(data=subset(data6, voteshare < 0.5 & fog2 == 0),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.man$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
  + geom_smooth(data=subset(data6, voteshare < 0.5 & fog2 == 1),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.may$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
  + geom_point(data=bin2.df, aes(x=mid, y=bin.mean, size=n), pch=1)
  + geom_point(data=bin2.dfb, aes(x=mid, y=bin.mean, size=n), pch=1,col="red")
  + geom_vline(xintercept=0)
  + geom_hline(yintercept=0, lty='dotted')
  + coord_cartesian(ylim = c(-0.05,1.05))
  + scale_x_continuous(breaks=seq(0,1,0.01),
                       limits=c(0.5-cct.run.win.may$h,
                                0.5+cct.run.win.may$h))
  
  + labs(size = '# Obs. in\n2-unit Bin')
  + theme_bw()
  + labs(x = "Vote share time t",
         y = "Prob. Running + Winning, time t+1")
  #  + ggtitle("Running + Winning by Form of Government")
  + annotate('text', x = 0.455, y = 0.28, 0.05, hjust=0, parse=TRUE, col="black",
             label = paste('hat(tau)==',
                           round(cct.run.win.man$coef['Conventional', ], 2),
                           '~(list(', round(cct.run.win.man$ci['Robust', 1], 2),
                           ',', round(cct.run.win.man$ci['Robust', 2], 2),
                           '))'))
  + annotate('text', x = 0.455, y = 0.6, 0.05, hjust=0, parse=TRUE, col="red",
             label = paste('hat(tau)==',
                           round(cct.run.win.may$coef['Conventional', ], 2),
                           '~(list(', round(cct.run.win.may$ci['Robust', 1], 2),
                           ',', round(cct.run.win.may$ci['Robust', 2], 2),
                           '))'))
  + annotate('text', x = 0.455, y = 0.33, 0.05, hjust=0, parse=TRUE,
             label = "Council-Manager",col="black")
  + annotate('text', x = 0.455, y = 0.65, 0.05, hjust=0, parse=TRUE,
             label = "Strong~Mayor",col="red")
  
)
dev.off()
######################



## By Term Limits:
data6$termlimits <- data$termlimits[match(data6$elec_index_next,data$elec_index)]
table(data6$termlimits) # worked
summarize(data,
          limits = length(unique(data$fips[data$termlimits==1])),
          nolimits = length(unique(data$fips[data$termlimits==0])))
# term length:
summarize(data,
          lengthmean = mean(data$termlength,na.rm=T)) # 2.9352
data6$termlength <- data$termlength[match(data6$elec_index_next,data$elec_index)]
length(unique(data6$elec_index[data6$termlength %in% c(1,3)])) # 433 elections
length(unique(data6$fips[which(data6$termlength ==1 | data6$termlength==3)])) # 58 cities
length(unique(data6$fips[which(data6$termlength ==1)])) # 18 cities
length(unique(data6$fips[which(data6$termlength ==2)])) # 387 cities
length(unique(data6$fips[which(data6$termlength ==3)])) # 41 cities
length(unique(data6$fips[which(data6$termlength ==4)])) # 542 cities


cct.run.termlimits <- with(data6[data6$termlimits==1,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                                 bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.notermlimits <- with(data6[data6$termlimits==0,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                                   bwselect="CCT", all=TRUE, kernel='uni'))

(Z_termlimits = (cct.run.notermlimits$coef[1,1] - cct.run.termlimits$coef[1,1])/sqrt(cct.run.notermlimits$se[3,1]^2 + cct.run.termlimits$se[3,1]^2))
SE_termlimits <- sqrt(cct.run.notermlimits$se[1,1]^2 + cct.run.termlimits$se[1,1]^2)
(Pvalue_termlimits <- 2*pnorm(-abs(Z_termlimits))) #

cct.run.win.termlimits <- with(data6[data6$termlimits==1,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                     bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.notermlimits <- with(data6[data6$termlimits==0,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                       bwselect="CCT", all=TRUE, kernel='uni'))

#####################
##### Table D-9: ####
#####################
rd.export5(list(cct.run.termlimits,cct.run.notermlimits,cct.run.win.termlimits,cct.run.win.notermlimits))
######################

(Z_termlimits = (cct.run.win.notermlimits$coef[1,1] - cct.run.win.termlimits$coef[1,1])/sqrt(cct.run.win.notermlimits$se[3,1]^2 + cct.run.win.termlimits$se[3,1]^2))
SE_termlimits <- sqrt(cct.run.win.notermlimits$se[1,1]^2 + cct.run.win.termlimits$se[1,1]^2)
(Pvalue_termlimits <- 2*pnorm(-abs(Z_termlimits))) # 


######################
##### Figure D-6 #####
######################
### Term Limits:
width <- .005
data6 <- mutate(data6, bin=cut(voteshare, breaks=seq(0, 1, width)))
bin2.df <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.notermlimits$h) & data6$voteshare<=(0.5+cct.run.win.notermlimits$h) & data6$termlimits==0)],
                                      data6$bin[which(data6$voteshare>=(0.5-cct.run.win.notermlimits$h) & data6$voteshare<=(0.5+cct.run.win.notermlimits$h) & data6$termlimits==0)], mean, na.rm=TRUE),
                      mid=seq(0 + width/2, 1 - width/2, width),
                      n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.notermlimits$h) & data6$voteshare<=(0.5+cct.run.win.notermlimits$h) & data6$termlimits==0)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.notermlimits$h) & data6$voteshare<=(0.5+cct.run.win.notermlimits$h) & data6$termlimits==0)],
                               length))
bin2.dfb <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.termlimits$h) & data6$voteshare<=(0.5+cct.run.win.termlimits$h) & data6$termlimits==1)],
                                       data6$bin[which(data6$voteshare>=(0.5-cct.run.win.termlimits$h) & data6$voteshare<=(0.5+cct.run.win.termlimits$h) & data6$termlimits==1)], mean, na.rm=TRUE),
                       mid=seq(0 + width/2, 1 - width/2, width),
                       n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.termlimits$h) & data6$voteshare<=(0.5+cct.run.win.termlimits$h) & data6$termlimits==1)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.termlimits$h) & data6$voteshare<=(0.5+cct.run.win.termlimits$h) & data6$termlimits==1)],
                                length))
tri <- function (x, h, c=0.5) pmax(0, 1 - abs((x - c) / h))

pdf("mayoral_rd_runwin_termlimits.pdf", width=8, height=6)
(ggplot(data6)
  + geom_point(data=subset(data6,termlimits == 0),
               aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5)
  + geom_point(data=subset(data6,termlimits == 1),
               aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5,col="red")
  + geom_smooth(data=subset(data6, voteshare > 0.5 & termlimits == 0),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.notermlimits$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
  + geom_smooth(data=subset(data6, voteshare > 0.5 & termlimits == 1),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.termlimits$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
  + geom_smooth(data=subset(data6, voteshare < 0.5 & termlimits == 0),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.notermlimits$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
  + geom_smooth(data=subset(data6, voteshare < 0.5 & termlimits == 1),
                aes(x = voteshare, y = run_and_win_next,
                    weight=tri(voteshare, cct.run.win.termlimits$h)),
                method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
  + geom_point(data=bin2.df, aes(x=mid, y=bin.mean, size=n), pch=1)
  + geom_point(data=bin2.dfb, aes(x=mid, y=bin.mean, size=n), pch=1,col="red")
  + geom_vline(xintercept=0)
  + geom_hline(yintercept=0, lty='dotted')
  + coord_cartesian(ylim = c(-0.05,1.05))
  + scale_x_continuous(breaks=seq(0,1,0.01),
                       limits=c(0.5-cct.run.win.termlimits$h,
                                0.5+cct.run.win.termlimits$h))
  
  + labs(size = '# Obs. in\n2-unit Bin')
  + theme_bw()
  + labs(x = "Vote share time t",
         y = "Prob. Running + Winning, time t+1")
  #  + ggtitle("Running + Winning by Term Limits")
  + annotate('text', x = 0.45, y = 0.28, 0.05, hjust=0, parse=TRUE, col="black",
             label = paste('hat(tau)==',
                           round(cct.run.win.notermlimits$coef['Conventional', ], 2),
                           '~(list(', round(cct.run.win.notermlimits$ci['Robust', 1], 2),
                           ',', round(cct.run.win.notermlimits$ci['Robust', 2], 2),
                           '))'))
  + annotate('text', x = 0.45, y = 0.6, 0.05, hjust=0, parse=TRUE, col="red",
             label = paste('hat(tau)==',
                           round(cct.run.win.termlimits$coef['Conventional', ], 2),
                           '~(list(', round(cct.run.win.termlimits$ci['Robust', 1], 2),
                           ',', round(cct.run.win.termlimits$ci['Robust', 2], 2),
                           '))'))
  + annotate('text', x = 0.45, y = 0.33, 0.05, hjust=0, parse=TRUE,
             label = "No~Term~Limits",col="black")
  + annotate('text', x = 0.45, y = 0.65, 0.05, hjust=0, parse=TRUE,
             label = "Term~Limits",col="red")
  
)
dev.off()
######################


#########
## Timing by other institutions:
########

# descriptives for switchers:

data$fog2 <- car::recode(data$fog,"1=1;2=0;else=NA") # managers = 0, strong mayor = 1?


sw.desc.tab <- data.frame(off=
                            rbind(
                                  round(mean(data$fog2[which(data$switcher==1 & data$concurrent==0)],na.rm=T),2), # 0.302
                                  paste0("(",round(sd(data$fog2[which(data$switcher==1 & data$concurrent==0)],na.rm=T),2),")"),
                                  round(mean(data$partisan[which(data$switcher==1 & data$concurrent==0)],na.rm=T),2), # 0.200
                                  paste0("(",round(sd(data$partisan[which(data$switcher==1 & data$concurrent==0)],na.rm=T),2),")"),
                                  round(mean(data$initiative[which(data$switcher==1 & data$concurrent==0)],na.rm=T),2), # 0.768
                                  paste0("(",round(sd(data$initiative[which(data$switcher==1 & data$concurrent==0)],na.rm=T),2),")"),
                                  round(mean(data$referendum[which(data$switcher==1 & data$concurrent==0)],na.rm=T),2), # 0.849
                                  paste0("(",round(sd(data$referendum[which(data$switcher==1 & data$concurrent==0)],na.rm=T),2),")"),
                                  length(data$fog2[which(data$switcher==1 & data$concurrent==0)])
                            ),
                          on=rbind(
                                   round(mean(data$fog2[which(data$switcher==1 & data$concurrent==1)],na.rm=T),2), # 0.353
                                   paste0("(",round(sd(data$fog2[which(data$switcher==1 & data$concurrent==1)],na.rm=T),2),")"),
                                   round(mean(data$partisan[which(data$switcher==1 & data$concurrent==1)],na.rm=T),2), # 0.257
                                   paste0("(",round(sd(data$partisan[which(data$switcher==1 & data$concurrent==1)],na.rm=T),2),")"),
                                   round(mean(data$initiative[which(data$switcher==1 & data$concurrent==1)],na.rm=T),2), # 0.651
                                   paste0("(",round(sd(data$initiative[which(data$switcher==1 & data$concurrent==1)],na.rm=T),2),")"),
                                   round(mean(data$referendum[which(data$switcher==1 & data$concurrent==1)],na.rm=T),2), # 0.879
                                   paste0("(",round(sd(data$referendum[which(data$switcher==1 & data$concurrent==1)],na.rm=T),2),")"),
                                   length(data$fog2[which(data$switcher==1 & data$concurrent==1)])
                          )
)


################
## Table D-10 ##
xtable(sw.desc.tab)
################


## Timing diff by partisan ballots:
cct.run.win.p_off <- with(data6[which(data6$partisan==1 & data6$concurrent_next==0),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.np_off <- with(data6[which(data6$partisan==0 & data6$concurrent_next==0),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.p_on <- with(data6[which(data6$partisan==1 & data6$concurrent_next==1),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.np_on <- with(data6[which(data6$partisan==0 & data6$concurrent_next==1),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,bwselect="CCT", all=TRUE, kernel='uni'))

######################
##### Table D-11: ####
######################
rd.export5(list(cct.run.win.p_on,cct.run.win.p_off,cct.run.win.np_on,cct.run.win.np_off))
######################

# significance test:
Z_timing_np = (cct.run.win.np_on$coef[1,1] - cct.run.win.np_off$coef[1,1])/sqrt(cct.run.win.np_on$se[3,1]^2 + cct.run.win.np_off$se[3,1]^2)
SE_timing_np <- sqrt(cct.run.win.np_on$se[1,1]^2 + cct.run.win.np_off$se[1,1]^2)
(Pvalue_timing_np <- 2*pnorm(-abs(Z_timing_np))) # 

Z_timing_p = (cct.run.win.p_on$coef[1,1] - cct.run.win.p_off$coef[1,1])/sqrt(cct.run.win.p_on$se[3,1]^2 + cct.run.win.p_off$se[3,1]^2)
SE_timing_p <- sqrt(cct.run.win.p_on$se[1,1]^2 + cct.run.win.p_off$se[1,1]^2)
(Pvalue_timing_p <- 2*pnorm(-abs(Z_timing_p))) # 



### timing differences by FOG:
cct.run.win.may_on <- with(data6[which(data6$fog==1 & data6$concurrent_next==1),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5))

cct.run.win.may_off <- with(data6[which(data6$fog==1 & data6$concurrent_next==0),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5))

cct.run.win.man_on <- with(data6[which(data6$fog==2 & data6$concurrent_next==1),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5))

cct.run.win.man_off <- with(data6[which(data6$fog==2 & data6$concurrent_next==0),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5))

(Z_timing_may = (cct.run.win.may_on$coef[1,1] - cct.run.win.may_off$coef[1,1])/sqrt(cct.run.win.may_on$se[3,1]^2 + cct.run.win.may_off$se[3,1]^2))
SE_timing_may <- sqrt(cct.run.win.may_on$se[1,1]^2 + cct.run.win.may_off$se[1,1]^2)
(Pvalue_timing_may <- 2*pnorm(-abs(Z_timing_may))) # 

(Z_timing_man = (cct.run.win.man_on$coef[1,1] - cct.run.win.man_off$coef[1,1])/sqrt(cct.run.win.man_on$se[3,1]^2 + cct.run.win.man_off$se[3,1]^2))
SE_timing_man <- sqrt(cct.run.win.man_on$se[1,1]^2 + cct.run.win.man_off$se[1,1]^2)
(Pvalue_timing_man <- 2*pnorm(-abs(Z_timing_man))) # 

######################
##### Table D-12: ####
######################
rd.export5(list(cct.run.win.may_on,cct.run.win.may_off,cct.run.win.man_on,cct.run.win.man_off))
######################


### Timing difference by term limits:
cct.run.notermlimits_off <- with(data6[which(data6$termlimits==0 & data6$concurrent_next==0),], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                                                                         bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.notermlimits_on <- with(data6[which(data6$termlimits==0 & data6$concurrent_next==1),], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                                                                        bwselect="CCT", all=TRUE, kernel='uni'))

(cct.run.notermlimits_off$coef[1,1] - cct.run.notermlimits_on$coef[1,1])
Z_timing_notermlimits2 = (cct.run.notermlimits_off$coef[1,1] - cct.run.notermlimits_on$coef[1,1])/sqrt(cct.run.notermlimits_off$se[3,1]^2 + cct.run.notermlimits_on$se[3,1]^2)
SE_timing_notermlimits2 <- sqrt(cct.run.notermlimits_off$se[1,1]^2 + cct.run.notermlimits_on$se[1,1]^2)
(Pvalue_timing_notermlimits2 <- 2*pnorm(-abs(Z_timing_notermlimits2))) # 


cct.run.win.notermlimits_off <- with(data6[which(data6$termlimits==0 & data6$concurrent_next==0),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                                                             bwselect="CCT", all=TRUE, kernel='uni'))

cct.run.win.notermlimits_on <- with(data6[which(data6$termlimits==0 & data6$concurrent_next==1),], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                                                                            bwselect="CCT", all=TRUE, kernel='uni'))

######################
##### Table D-13: ####
######################
rd.export5(list(cct.run.notermlimits_on,cct.run.notermlimits_off,cct.run.win.notermlimits_on,cct.run.win.notermlimits_off))
######################
(cct.run.win.notermlimits_off$coef[1,1] - cct.run.win.notermlimits_on$coef[1,1])
Z_timing_notermlimits = (cct.run.win.notermlimits_off$coef[1,1] - cct.run.win.notermlimits_on$coef[1,1])/sqrt(cct.run.win.notermlimits_off$se[3,1]^2 + cct.run.win.notermlimits_on$se[3,1]^2)
SE_timing_notermlimits <- sqrt(cct.run.win.notermlimits_off$se[1,1]^2 + cct.run.win.notermlimits_on$se[1,1]^2)
(Pvalue_timing_notermlimits <- 2*pnorm(-abs(Z_timing_notermlimits))) # 







######################
##### Appendix E #####
######################
# Impact of Media Environment:

### Incumbency advantage by presence of daily paper:
cct.run.paper <- with(data6[data6$dailypaper==1,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                   bwselect="CCT", all=TRUE, kernel='uni'))
cct.run.paper2 <- with(data6[data6$dailypaper==1,], rdrobust(y=run_next , x=voteshare , c=0.5))

cct.run.nopaper <- with(data6[data6$dailypaper==0,], rdrobust(y=run_next , x=voteshare , c=0.5,
                                                   bwselect="CCT", all=TRUE, kernel='uni'))
cct.run.nopaper2 <- with(data6[data6$dailypaper==0,], rdrobust(y=run_next , x=voteshare , c=0.5))


cct.run.win.paper <- with(data6[data6$dailypaper==1,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5))
                          
cct.run.win.nopaper <- with(data6[data6$dailypaper==0,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5,
                                                       bwselect="CCT", all=TRUE, kernel='uni'))
cct.run.win.nopaper2 <- with(data6[data6$dailypaper==0,], rdrobust(y=run_and_win_next , x=voteshare , c=0.5))

# significance test:
(cct.run.win.nopaper$coef[1,1] - cct.run.win.paper$coef[1,1])
Z_paper = (cct.run.win.nopaper$coef[1,1] - cct.run.win.paper$coef[1,1])/sqrt(cct.run.win.nopaper$se[3,1]^2 + cct.run.win.paper$se[3,1]^2)
SE_paper <- sqrt(cct.run.win.nopaper$se[1,1]^2 + cct.run.win.paper$se[1,1]^2)
(Pvalue_paper <- 2*pnorm(-abs(Z_paper))) # 

#######################
##### Table E-14: #####
#######################
rd.export5(list(cct.run.paper,cct.run.nopaper,cct.run.win.paper,cct.run.win.nopaper))
#######################



######################
##### Figure E-7 #####
######################
### Media Environment:
width <- .005
data6 <- mutate(data6, bin=cut(voteshare, breaks=seq(0, 1, width)))
bin2.df <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.nopaper$h) & data6$voteshare<=(0.5+cct.run.win.nopaper$h) & data6$dailypaper==0)],
                                      data6$bin[which(data6$voteshare>=(0.5-cct.run.win.nopaper$h) & data6$voteshare<=(0.5+cct.run.win.nopaper$h) & data6$dailypaper==0)], mean, na.rm=TRUE),
                      mid=seq(0 + width/2, 1 - width/2, width),
                      n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.nopaper$h) & data6$voteshare<=(0.5+cct.run.win.nopaper$h) & data6$dailypaper==0)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.nopaper$h) & data6$voteshare<=(0.5+cct.run.win.nopaper$h) & data6$dailypaper==0)],
                               length))
bin2.dfb <- data.frame(bin.mean=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.paper$h) & data6$voteshare<=(0.5+cct.run.win.paper$h) & data6$dailypaper==1)],
                                       data6$bin[which(data6$voteshare>=(0.5-cct.run.win.paper$h) & data6$voteshare<=(0.5+cct.run.win.paper$h) & data6$dailypaper==1)], mean, na.rm=TRUE),
                       mid=seq(0 + width/2, 1 - width/2, width),
                       n=tapply(data6$run_and_win_next[which(data6$voteshare>=(0.5-cct.run.win.paper$h) & data6$voteshare<=(0.5+cct.run.win.paper$h) & data6$dailypaper==1)], data6$bin[which(data6$voteshare>=(0.5-cct.run.win.paper$h) & data6$voteshare<=(0.5+cct.run.win.paper$h) & data6$dailypaper==1)],
                                length))
tri <- function (x, h, c=0.5) pmax(0, 1 - abs((x - c) / h))

pdf("mayoral_rd_runwin_media.pdf", width=8, height=6)
(ggplot(data6)
 + geom_point(data=subset(data6,dailypaper == 0),
              aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5)
 + geom_point(data=subset(data6,dailypaper == 1),
              aes(x = voteshare, y = run_and_win_next), shape="|", alpha=.5,col="red")
 + geom_smooth(data=subset(data6, voteshare > 0.5 & dailypaper == 0),
               aes(x = voteshare, y = run_and_win_next,
                   weight=tri(voteshare, cct.run.win.nopaper$h)),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
 + geom_smooth(data=subset(data6, voteshare > 0.5 & dailypaper == 1),
               aes(x = voteshare, y = run_and_win_next,
                   weight=tri(voteshare, cct.run.win.paper$h)),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
 + geom_smooth(data=subset(data6, voteshare < 0.5 & dailypaper == 0),
               aes(x = voteshare, y = run_and_win_next,
                   weight=tri(voteshare, cct.run.win.nopaper$h)),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="black")
 + geom_smooth(data=subset(data6, voteshare < 0.5 & dailypaper == 1),
               aes(x = voteshare, y = run_and_win_next,
                   weight=tri(voteshare, cct.run.win.paper$h)),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5, col="red")
 + geom_point(data=bin2.df, aes(x=mid, y=bin.mean, size=n), pch=1)
 + geom_point(data=bin2.dfb, aes(x=mid, y=bin.mean, size=n), pch=1,col="red")
 + geom_vline(xintercept=0)
 + geom_hline(yintercept=0, lty='dotted')
 + coord_cartesian(ylim = c(-0.05,1.05))
 + scale_x_continuous(breaks=seq(0,1,0.01),
                      limits=c(0.5-cct.run.win.paper$h,
                               0.5+cct.run.win.paper$h))
 
 + labs(size = '# Obs. in\n2-unit Bin')
 + theme_bw()
 + labs(x = "Vote share time t",
        y = "Prob. Running + Winning, time t+1")
#  + ggtitle("Running + Winning by Media Environment")
 + annotate('text', x = 0.455, y = 0.28, 0.05, hjust=0, parse=TRUE, col="black",
            label = paste('hat(tau)==',
                          round(cct.run.win.nopaper$coef['Conventional', ], 2),
                          '~(list(', round(cct.run.win.nopaper$ci['Robust', 1], 2),
                          ',', round(cct.run.win.nopaper$ci['Robust', 2], 2),
                          '))'))
 + annotate('text', x = 0.455, y = 0.6, 0.05, hjust=0, parse=TRUE, col="red",
            label = paste('hat(tau)==',
                          round(cct.run.win.paper$coef['Conventional', ], 2),
                          '~(list(', round(cct.run.win.paper$ci['Robust', 1], 2),
                          ',', round(cct.run.win.paper$ci['Robust', 2], 2),
                          '))'))
 + annotate('text', x = 0.455, y = 0.33, 0.05, hjust=0, parse=TRUE,
            label = "No~Daily~Paper",col="black")
 + annotate('text', x = 0.455, y = 0.65, 0.05, hjust=0, parse=TRUE,
            label = "Daily~Paper",col="red")
 
)
dev.off()
######################



#######################
##### Appendix F: #####
#######################



### continuity test comparing before/after treatment changing election timing:
bw <- cct.run$h

check_data_post <- data6[which((data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta2==1 & data6$concurrent_next==1) | (data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==1 & data6$concurrent_next==1)),]

width <- .005
check_data_post <- mutate(check_data_post, bin=cut(voteshare_adj, breaks=seq(-0.5, 0.5, width)))
bin2.df.DID.post <- data.frame(bin.mean=tapply(check_data_post$voteshare_adj,
                                               check_data_post$bin, mean, na.rm=TRUE),
                               mid=seq(-0.5 + width/2, 0.5 - width/2, width),
                               n=tapply(check_data_post$voteshare_adj, check_data_post$bin,
                                        length))
mc.rd.DID.post <- lm(n ~ mid * (mid>=0), data=bin2.df.DID.post[which(bin2.df.DID.post$mid>=(-bw) & bin2.df.DID.post$mid<=(bw)),])
summary(mc.rd.DID.post) # beta = 0.015, p=0.99

check_data_pre <- data6[which((data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead==1 & data6$concurrent_next==0) | (data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead2==1 & data6$concurrent_delta_lead==0 & data6$concurrent_next==0)),]

check_data_pre <- mutate(check_data_pre, bin=cut(voteshare_adj, breaks=seq(-0.5, 0.5, width)))
bin2.df.DID.pre <- data.frame(bin.mean=tapply(check_data_pre$voteshare_adj,
                                              check_data_pre$bin, mean, na.rm=TRUE),
                              mid=seq(-0.5 + width/2, 0.5 - width/2, width),
                              n=tapply(check_data_pre$voteshare_adj, check_data_pre$bin,
                                       length))
mc.rd.DID.pre <- lm(n ~ mid * (mid>=0), data=bin2.df.DID.pre[which(bin2.df.DID.pre$mid>=(-bw) & bin2.df.DID.pre$mid<=(bw)),])
summary(mc.rd.DID.pre) # beta = 0.000, p=1

bin2.df.DID.diff <- data.frame(bin.mean=tapply(check_data_pre$voteshare_adj,
                                               check_data_pre$bin, mean, na.rm=TRUE),
                               mid=seq(-0.5 + width/2, 0.5 - width/2, width),
                               n=(bin2.df.DID.post$n-bin2.df.DID.pre$n)
)
mc.rd.DID.diff <- lm(n ~ mid * (mid>=0), data=bin2.df.DID.diff[which(bin2.df.DID.diff$mid>=(-bw) & bin2.df.DID.diff$mid<=(bw)),])
summary(mc.rd.DID.diff) # beta = 0.015, p=0.99

######################
##### Table F-15 #####
######################
xtable(rbind(summary(mc.rd.DID.diff)$coefficients[3,],summary(mc.rd.DID.pre)$coefficients[3,],summary(mc.rd.DID.post)$coefficients[3,]))
######################


## Diff-in-Disc:

# set bandwidth
bw <- cct.run.win.on$h


fit_postlose <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==-1),])) # 40 obs
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==-1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
postlose_se <- sd(allcoefs) # SD
postlose_lo <- quantile(allcoefs,0.025)
postlose_hi <- quantile(allcoefs,0.975)

fit_post2lose <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta2==-1),])) # only 38 obs
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta2==-1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
post2lose_se <- sd(allcoefs) # SD
post2lose_lo <- quantile(allcoefs,0.025)
post2lose_hi <- quantile(allcoefs,0.975)

fit_prelose <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead==-1 & data6$concurrent_delta==0 & data6$concurrent_next==1),])) # 28 obs
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead==-1 & data6$concurrent_delta==0 & data6$concurrent_next==1),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
prelose_se <- sd(allcoefs) # SD
prelose_lo <- quantile(allcoefs,0.025)
prelose_hi <- quantile(allcoefs,0.975)

fit_pre2lose <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead2==-1 & data6$concurrent_delta_lead==0 & data6$concurrent_delta==0 & data6$concurrent_next==1),])) # 18 obs
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead2==-1 & data6$concurrent_delta_lead==0 & data6$concurrent_delta==0 & data6$concurrent_next==1),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
pre2lose_se <- sd(allcoefs) # SD
pre2lose_lo <- quantile(allcoefs,0.025)
pre2lose_hi <- quantile(allcoefs,0.975)


inc_df <- data.frame(val=c(fit_pre2lose$coef[3],fit_prelose$coef[3],fit_postlose$coef[3],fit_post2lose$coef[3],fit_pre2on$coef[3],fit_preon$coef[3],fit_poston$coef[3],fit_post2on$coef[3],fit_pre2off$coef[3],fit_preoff$coef[3],fit_postoff$coef[3],fit_post2off$coef[3]),time=c("2 periods pre-switch","pre-switch","post-switch","2 periods post-switch","2 periods pre-switch","pre-switch","post-switch","2 periods post-switch","2 periods pre-switch","pre-switch","post-switch","2 periods post-switch"),group=c("move to off-cycle","move to off-cycle","move to off-cycle","move to off-cycle","stay on-cycle","stay on-cycle","stay on-cycle","stay on-cycle","stay off-cycle","stay off-cycle","stay off-cycle","stay off-cycle"))
inc_df <- data.frame(val=c(fit_pre2lose$coef[3],fit_prelose$coef[3],fit_postlose$coef[3],fit_post2lose$coef[3],fit_pre2on$coef[3],fit_preon$coef[3],fit_poston$coef[3],fit_post2on$coef[3],fit_pre2off$coef[3],fit_preoff$coef[3],fit_postoff$coef[3],fit_post2off$coef[3]),se=c(pre2lose_se,prelose_se,postlose_se,post2lose_se,pre2on_se,preon_se,poston_se,post2on_se,pre2off_se,preoff_se,postoff_se,post2off_se),ci.lo=c(pre2lose_lo,prelose_lo,postlose_lo,post2lose_lo,pre2on_lo,preon_lo,poston_lo,post2on_lo,pre2off_lo,preoff_lo,postoff_lo,post2off_lo),ci.hi=c(pre2lose_hi,prelose_hi,postlose_hi,post2lose_hi,pre2on_hi,preon_hi,poston_hi,post2on_hi,pre2off_hi,preoff_hi,postoff_hi,post2off_hi),time=c("2 periods pre-switch","pre-switch","post-switch","2 periods post-switch","2 periods pre-switch","pre-switch","post-switch","2 periods post-switch","2 periods pre-switch","pre-switch","post-switch","2 periods post-switch"),group=c("move to off-cycle","move to off-cycle","move to off-cycle","move to off-cycle","stay on-cycle","stay on-cycle","stay on-cycle","stay on-cycle","stay off-cycle","stay off-cycle","stay off-cycle","stay off-cycle"))

# inc_df$ci.lo <- inc_df$val+qnorm(0.025)*inc_df$se
# inc_df$ci.hi <- inc_df$val+qnorm(0.975)*inc_df$se
timeorder <- c("2 periods pre-switch","pre-switch","post-switch","2 periods post-switch")
pdf("diff_in_disc_timingLOSE_noCI.pdf",height=6,width=8)
ggplot(inc_df,aes(x=time,y=val,col=group,shape=group)) +
  #   geom_errorbar(aes(x = time, y = val, ymin = ci.lo, ymax = ci.hi,col=group), size = 1.5, width = 0.2,alpha=0.5) +
  geom_smooth(method="loess",aes(group=group)) +
  geom_point(size=4) +
  scale_y_continuous("Incumbency Advantage (Pr. Running + Winning)") + scale_x_discrete("Time",limits=timeorder) +
  coord_cartesian(ylim=c(-0.1,1.1)) + 
  ggtitle("Incumbency Advantage,\nChanges to Election Timing")
dev.off()



## what about cities that move TO concurrent?
fit_postgain <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==1 & data6$concurrent_next==1),]))
# fit_postgain_out <- coeftest(fit_postgain,vcov=vcovCluster(fit_postgain,cluster=data6$gov_id[which(data6$voteshare_adj>=-0.05 & data6$voteshare_adj<=0.05 & data6$concurrent_delta==1 & !is.na(data6$run_and_win_next))])) #  significant still, but only 58 clusters... should block-bootstrap
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
postgain_se <- sd(allcoefs) # SD
postgain_lo <- quantile(allcoefs,0.025)
postgain_hi <- quantile(allcoefs,0.975)
fit_postgain$coef[3]+c(qnorm(0.025),qnorm(0.975))*postgain_se # CIs


fit_post2gain <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta2==1 & data6$concurrent_next==1),])) # only 46 obs

# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta2==1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
post2gain_se <- sd(allcoefs) # SD
post2gain_lo <- quantile(allcoefs,0.025)
post2gain_hi <- quantile(allcoefs,0.975)
fit_post2gain$coef[3]+c(qnorm(0.025),qnorm(0.975))*post2gain_se # CIs


fit_pregain <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead==1 & data6$concurrent_next==0),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead==1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
pregain_se <- sd(allcoefs) # SD
pregain_lo <- quantile(allcoefs,0.025)
pregain_hi <- quantile(allcoefs,0.975)
fit_pregain$coef[3]+c(qnorm(0.025),qnorm(0.975))*pregain_se # CIs

fit_pre2gain <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead2==1 & data6$concurrent_delta_lead==0 & data6$concurrent_next==0),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead2==1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
pre2gain_se <- sd(allcoefs) # SD
pre2gain_lo <- quantile(allcoefs,0.025)
pre2gain_hi <- quantile(allcoefs,0.975)
fit_pre2gain$coef[3]+c(qnorm(0.025),qnorm(0.975))*pre2gain_se # CIs

fit_poston <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==0 & data6$concurrent_next==1),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==0 & data6$concurrent_next==1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
poston_se <- sd(allcoefs) # SD
poston_lo <- quantile(allcoefs,0.025)
poston_hi <- quantile(allcoefs,0.975)
fit_poston$coef[3]+c(qnorm(0.025),qnorm(0.975))*poston_se # CIs

fit_post2on <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta2==0 & data6$concurrent_next==1),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta2==0 & data6$concurrent_next==1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
post2on_se <- sd(allcoefs) # SD
post2on_lo <- quantile(allcoefs,0.025)
post2on_hi <- quantile(allcoefs,0.975)
fit_post2on$coef[3]+c(qnorm(0.025),qnorm(0.975))*post2on_se # CIs

fit_preon <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==0 & data6$concurrent_delta_lead==0 & data6$concurrent_next==1),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead==0 & data6$concurrent_next==1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
preon_se <- sd(allcoefs) # SD
preon_lo <- quantile(allcoefs,0.025)
preon_hi <- quantile(allcoefs,0.975)
fit_preon$coef[3]+c(qnorm(0.025),qnorm(0.975))*preon_se # CIs

fit_pre2on <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==0 & data6$concurrent_delta2==0 & data6$concurrent_delta_lead==0 & data6$concurrent_delta_lead2==0 & data6$concurrent_next==1),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead2==0 & data6$concurrent_next==1 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
pre2on_se <- sd(allcoefs) # SD
pre2on_lo <- quantile(allcoefs,0.025)
pre2on_hi <- quantile(allcoefs,0.975)
fit_pre2on$coef[3]+c(qnorm(0.025),qnorm(0.975))*pre2on_se # CIs

fit_postoff <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==0 & data6$concurrent_next==0),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==0 & data6$concurrent_next==0 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
postoff_se <- sd(allcoefs) # SD
postoff_lo <- quantile(allcoefs,0.025)
postoff_hi <- quantile(allcoefs,0.975)
fit_postoff$coef[3]+c(qnorm(0.025),qnorm(0.975))*postoff_se # CIs

fit_post2off <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta2==0 & data6$concurrent_next==0),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta2==0 & data6$concurrent_next==0 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
post2off_se <- sd(allcoefs) # SD
post2off_lo <- quantile(allcoefs,0.025)
post2off_hi <- quantile(allcoefs,0.975)
fit_post2off$coef[3]+c(qnorm(0.025),qnorm(0.975))*post2off_se # CIs


fit_preoff <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==0 & data6$concurrent_delta_lead==0 & data6$concurrent_next==0),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead==0 & data6$concurrent_next==0 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
preoff_se <- sd(allcoefs) # SD
preoff_lo <- quantile(allcoefs,0.025)
preoff_hi <- quantile(allcoefs,0.975)
fit_preoff$coef[3]+c(qnorm(0.025),qnorm(0.975))*preoff_se # CIs

fit_pre2off <- (lm(run_and_win_next ~ voteshare_adj * (voteshare>=0.5),data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta==0 & data6$concurrent_delta_lead==0 & data6$concurrent_delta_lead2==0 & data6$concurrent_next==0),]))
# Block bootstrap:
d <- data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw & data6$concurrent_delta_lead2==0 & data6$concurrent_next==0 & !is.na(data6$run_and_win_next)),]
iter = 10000 #number of iterations
# mat = matrix(NA, nrow=nrow(d), ncol=iter)
allcoefs <- NULL
for(i in 1:iter){
  block.boot = sample(unique(d$gov_id),length(unique(d$gov_id)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$gov_id %in% block.boot[j]))
  }
  thisfit <- lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=d[boot,])
  allcoefs <- c(allcoefs,thisfit$coef[3])
}
pre2off_se <- sd(allcoefs) # SD
pre2off_lo <- quantile(allcoefs,0.025)
pre2off_hi <- quantile(allcoefs,0.975)
fit_pre2off$coef[3]+c(qnorm(0.025),qnorm(0.975))*pre2off_se # CIs

inc_df <- data.frame(val=c(fit_pre2gain$coef[3],fit_pregain$coef[3],fit_postgain$coef[3],fit_post2gain$coef[3],fit_pre2on$coef[3],fit_preon$coef[3],fit_poston$coef[3],fit_post2on$coef[3],fit_pre2off$coef[3],fit_preoff$coef[3],fit_postoff$coef[3],fit_post2off$coef[3]),time=c("2 periods pre-switch","pre-switch","post-switch","2 periods post-switch","2 periods pre-switch","pre-switch","post-switch","2 periods post-switch","2 periods pre-switch","pre-switch","post-switch","2 periods post-switch"),group=c("move to concurrent","move to concurrent","move to concurrent","move to concurrent","stay concurrent","stay concurrent","stay concurrent","stay concurrent","stay off-cycle","stay off-cycle","stay off-cycle","stay off-cycle"))
inc_df <- data.frame(val=c(fit_pre2gain$coef[3],fit_pregain$coef[3],fit_postgain$coef[3],fit_post2gain$coef[3],fit_pre2on$coef[3],fit_preon$coef[3],fit_poston$coef[3],fit_post2on$coef[3],fit_pre2off$coef[3],fit_preoff$coef[3],fit_postoff$coef[3],fit_post2off$coef[3]),se=c(pre2gain_se,pregain_se,postgain_se,post2gain_se,pre2on_se,preon_se,poston_se,post2on_se,pre2off_se,preoff_se,postoff_se,post2off_se),ci.lo=c(pre2gain_lo,pregain_lo,postgain_lo,post2gain_lo,pre2on_lo,preon_lo,poston_lo,post2on_lo,pre2off_lo,preoff_lo,postoff_lo,post2off_lo),ci.hi=c(pre2gain_hi,pregain_hi,postgain_hi,post2gain_hi,pre2on_hi,preon_hi,poston_hi,post2on_hi,pre2off_hi,preoff_hi,postoff_hi,post2off_hi),time=c("2 periods pre-switch","pre-switch","post-switch","2 periods post-switch","2 periods pre-switch","pre-switch","post-switch","2 periods post-switch","2 periods pre-switch","pre-switch","post-switch","2 periods post-switch"),group=c("move to concurrent","move to concurrent","move to concurrent","move to concurrent","stay concurrent","stay concurrent","stay concurrent","stay concurrent","stay off-cycle","stay off-cycle","stay off-cycle","stay off-cycle"))
# inc_df$ci.lo <- inc_df$val+qnorm(0.025)*inc_df$se
# inc_df$ci.hi <- inc_df$val+qnorm(0.975)*inc_df$se

timeorder <- c("2 periods pre-switch","pre-switch","post-switch","2 periods post-switch")


pdf("diff_in_disc_timing_noCI.pdf",height=6,width=8)
ggplot(inc_df,aes(x=time,y=val,col=group,shape=group)) +
#   geom_errorbar(aes(x = time, y = val, ymin = ci.lo, ymax = ci.hi,col=group), size = 1.5, width = 0.2,alpha=0.5) +
  geom_smooth(method="loess",aes(group=group)) +
  geom_point(size=4) +
  scale_y_continuous("Incumbency Advantage\n(Probability Running + Winning)") + scale_x_discrete("Time",limits=timeorder) +
  #   ggtitle("Incumbency Advantage,\nChanges to Election Timing") +
  coord_cartesian(ylim=c(-0.1,1.1))
dev.off()

## Table:
did_tab <- cbind((c(rep("$t-2$",3),rep("$t-1$",3),rep("$t$",3),rep("$t+1$",3))),c("Will become on-cycle","Will stay on-cycle","Will stay off-cycle","Will become on-cycle","Will stay on-cycle","Will stay off-cycle","Became on-cycle","Stayed on-cycle","Stayed off-cycle","Became on-cycle","Stayed on-cycle","Stayed off-cycle"),round(c(inc_df$val[inc_df$time=="2 periods pre-switch"],inc_df$val[inc_df$time=="pre-switch"],inc_df$val[inc_df$time=="post-switch"],inc_df$val[inc_df$time=="2 periods post-switch"]),2),c(paste0("(",round(inc_df$ci.lo[inc_df$time=="2 periods pre-switch"],2),", ",round(inc_df$ci.hi[inc_df$time=="2 periods pre-switch"],2),")"),paste0("(",round(inc_df$ci.lo[inc_df$time=="pre-switch"],2),", ",round(inc_df$ci.hi[inc_df$time=="pre-switch"],2),")"),paste0("(",round(inc_df$ci.lo[inc_df$time=="post-switch"],2),", ",round(inc_df$ci.hi[inc_df$time=="post-switch"],2),")"),paste0("(",round(inc_df$ci.lo[inc_df$time=="2 periods post-switch"],2),", ",round(inc_df$ci.hi[inc_df$time=="2 periods post-switch"],2),")")),c(nrow(model.frame(fit_pre2gain)),nrow(model.frame(fit_pre2on)),nrow(model.frame(fit_pre2off)),nrow(model.frame(fit_pregain)),nrow(model.frame(fit_preon)),nrow(model.frame(fit_preoff)),nrow(model.frame(fit_postgain)),nrow(model.frame(fit_poston)),nrow(model.frame(fit_postoff)),nrow(model.frame(fit_post2gain)),nrow(model.frame(fit_post2on)),nrow(model.frame(fit_post2off))),rep(round(cct.run.win.on$h,2),times=12))

print(xtable(did_tab),include.rownames=F)

## DID estimates:
# store estimates
did.on2 <- (fit_post2gain$coef[3]-fit_postgain$coef[3])-(fit_post2on$coef[3]-fit_poston$coef[3]) # DID estimate of 0.45
did.off2 <- (fit_post2gain$coef[3]-fit_postgain$coef[3])-(fit_post2off$coef[3]-fit_postoff$coef[3]) # DID estimate of 0.37
did.on <- (fit_postgain$coef[3]-fit_pregain$coef[3])-(fit_poston$coef[3]-fit_preon$coef[3]) # DID estimate of -0.111
did.off <- (fit_postgain$coef[3]-fit_pregain$coef[3])-(fit_postoff$coef[3]-fit_preoff$coef[3]) # DID estimate of -0.104
did.on.sep <- (fit_post2gain$coef[3]-fit_pregain$coef[3])-(fit_post2on$coef[3]-fit_preon$coef[3]) # DID estimate of 0.320
did.off.sep <- (fit_post2gain$coef[3]-fit_pregain$coef[3])-(fit_post2off$coef[3]-fit_preoff$coef[3]) # DID estimate of 0.315
did.lose.on.sep <- (fit_post2lose$coef[3]-fit_prelose$coef[3])-(fit_post2on$coef[3]-fit_preon$coef[3]) # DID estimate of -0.678
did.lose.off.sep <- (fit_post2lose$coef[3]-fit_prelose$coef[3])-(fit_post2off$coef[3]-fit_preoff$coef[3]) # DID estimate of -0.683
plac.on <- (fit_pregain$coef[3]-fit_pre2gain$coef[3])-(fit_preon$coef[3]-fit_pre2on$coef[3]) # placebo est of -0.005
plac.off <- (fit_pregain$coef[3]-fit_pre2gain$coef[3])-(fit_preoff$coef[3]-fit_pre2off$coef[3]) # 0.158

# Block bootstrap:
start <- Sys.time()
d <- data6
iter = 10000 #number of iterations
alldiffs.gain.on2 <- NULL
alldiffs.gain.off2 <- NULL
alldiffs.gain.on <- NULL
alldiffs.gain.off <- NULL
alldiffs.gain.on.sep <- NULL
alldiffs.gain.off.sep <- NULL
allplac.gain.on <- NULL
allplac.gain.off <- NULL
alldiffs.lose.on2 <- NULL
alldiffs.lose.off2 <- NULL
alldiffs.lose.on <- NULL
alldiffs.lose.off <- NULL
alldiffs.lose.on.sep <- NULL
alldiffs.lose.off.sep <- NULL
allplac.lose.on <- NULL
allplac.lose.off <- NULL
set.seed(02139)
for(i in 1:iter){
  block.boot = sample(unique(d$fips),length(unique(d$fips)),replace=T)
  boot = c()
  for(j in 1:length(block.boot)){
    boot = c(boot,which(d$fips %in% block.boot[j]))
  }
  # create bootstrapped data:
  thisdata <- d[boot,]
  
  # run regressions:
  fit_post2gain <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta2==1),]))
  fit_postgain <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta==1),]))
  fit_pregain <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta_lead==1),]))
  fit_pre2gain <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta_lead2==1),]))
  fit_post2lose <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta2==-1),]))
  fit_postlose <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta==-1),]))
  fit_prelose <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta_lead==-1),]))
  fit_pre2lose <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta_lead2==-1),]))
  fit_post2on <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta2==0 & thisdata$concurrent_next==1),]))
  fit_poston <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta==0 & thisdata$concurrent_next==1),]))
  fit_preon <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta_lead==0 & thisdata$concurrent_next==1),]))
  fit_pre2on <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta_lead2==0 & thisdata$concurrent_next==1),]))
  fit_post2off <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta2==0 & thisdata$concurrent_next==0),]))
  fit_postoff <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta==0 & thisdata$concurrent_next==0),]))
  fit_preoff <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta_lead==0 & thisdata$concurrent_next==0),]))
  fit_pre2off <- (lm(run_and_win_next ~ voteshare_adj * (voteshare_adj>=0),data=thisdata[which(thisdata$voteshare_adj>=-bw & thisdata$voteshare_adj<=bw & thisdata$concurrent_delta_lead2==0 & thisdata$concurrent_next==0),]))
  
  # store quantity of interest:
  alldiffs.gain.on2 <- c(alldiffs.gain.on2,(fit_post2gain$coef[3]-fit_postgain$coef[3])-(fit_post2on$coef[3]-fit_poston$coef[3]))
  alldiffs.gain.off2 <- c(alldiffs.gain.off2,(fit_post2gain$coef[3]-fit_postgain$coef[3])-(fit_post2off$coef[3]-fit_postoff$coef[3]))
  alldiffs.gain.on <- c(alldiffs.gain.on,(fit_postgain$coef[3]-fit_pregain$coef[3])-(fit_poston$coef[3]-fit_preon$coef[3]))
  alldiffs.gain.off <- c(alldiffs.gain.off,(fit_postgain$coef[3]-fit_pregain$coef[3])-(fit_postoff$coef[3]-fit_preoff$coef[3]))
  alldiffs.gain.on.sep <- c(alldiffs.gain.on.sep,(fit_post2gain$coef[3]-fit_pregain$coef[3])-(fit_post2on$coef[3]-fit_preon$coef[3]))
  alldiffs.gain.off.sep <- c(alldiffs.gain.off.sep,(fit_post2gain$coef[3]-fit_pregain$coef[3])-(fit_post2off$coef[3]-fit_preoff$coef[3]))
  allplac.gain.on <- c(allplac.gain.on,(fit_pregain$coef[3]-fit_pre2gain$coef[3])-(fit_preon$coef[3]-fit_pre2on$coef[3]))
  allplac.gain.off <- c(allplac.gain.off,(fit_pregain$coef[3]-fit_pre2gain$coef[3])-(fit_preoff$coef[3]-fit_pre2off$coef[3]))
  alldiffs.lose.on2 <- c(alldiffs.lose.on2,(fit_post2lose$coef[3]-fit_postlose$coef[3])-(fit_post2on$coef[3]-fit_poston$coef[3]))
  alldiffs.lose.off2 <- c(alldiffs.lose.off2,(fit_post2lose$coef[3]-fit_postlose$coef[3])-(fit_post2off$coef[3]-fit_postoff$coef[3]))
  alldiffs.lose.on <- c(alldiffs.lose.on,(fit_postlose$coef[3]-fit_prelose$coef[3])-(fit_poston$coef[3]-fit_preon$coef[3]))
  alldiffs.lose.off <- c(alldiffs.lose.off,(fit_postlose$coef[3]-fit_prelose$coef[3])-(fit_postoff$coef[3]-fit_preoff$coef[3]))
  alldiffs.lose.on.sep <- c(alldiffs.lose.on.sep,(fit_post2lose$coef[3]-fit_prelose$coef[3])-(fit_post2on$coef[3]-fit_preon$coef[3]))
  alldiffs.lose.off.sep <- c(alldiffs.lose.off.sep,(fit_post2lose$coef[3]-fit_prelose$coef[3])-(fit_post2off$coef[3]-fit_preoff$coef[3]))
  allplac.lose.on <- c(allplac.lose.on,(fit_prelose$coef[3]-fit_pre2lose$coef[3])-(fit_preon$coef[3]-fit_pre2on$coef[3]))
  allplac.lose.off <- c(allplac.lose.off,(fit_prelose$coef[3]-fit_pre2lose$coef[3])-(fit_preoff$coef[3]-fit_pre2off$coef[3]))
}
end <- Sys.time()
end-start # took 1.86 hours
sd(alldiffs.gain.on)
sd(alldiffs.gain.off)
ci.did.gain.on2 <- quantile(alldiffs.gain.on2,c(0.025,0.975),na.rm=T)
ci.did.gain.off2 <- quantile(alldiffs.gain.off2,c(0.025,0.975),na.rm=T)
# ci.did.gain.on2.10 <- quantile(alldiffs.gain.on2,c(0.05,0.95),na.rm=T)
ci.did.gain.off2.10 <- quantile(alldiffs.gain.off2,c(0.05,0.95),na.rm=T)

ci.did.gain.on <- quantile(alldiffs.gain.on,c(0.025,0.975),na.rm=T)
ci.did.gain.off <- quantile(alldiffs.gain.off,c(0.025,0.975),na.rm=T)
ci.did.gain.on.sep <- quantile(alldiffs.gain.on.sep,c(0.025,0.975),na.rm=T)
ci.did.gain.off.sep <- quantile(alldiffs.gain.off.sep,c(0.025,0.975),na.rm=T)

ci.did.lose.on.sep <- quantile(alldiffs.lose.on.sep,c(0.025,0.975),na.rm=T)
ci.did.lose.off.sep <- quantile(alldiffs.lose.off.sep,c(0.025,0.975),na.rm=T)

ci.plac.gain.on <- quantile(allplac.gain.on,c(0.025,0.975),na.rm=T)
ci.plac.gain.off <- quantile(allplac.gain.off,c(0.025,0.975),na.rm=T)
ci.plac.lose.on <- quantile(allplac.lose.on,c(0.025,0.975),na.rm=T)
ci.plac.lose.off <- quantile(allplac.lose.off,c(0.025,0.975),na.rm=T)
# at 90% levels:
ci.plac.gain.on90 <- quantile(allplac.gain.on,c(0.05,0.95),na.rm=T)
ci.plac.gain.off90 <- quantile(allplac.gain.off,c(0.05,0.95),na.rm=T)
ci.plac.lose.on90 <- quantile(allplac.lose.on,c(0.05,0.95),na.rm=T)
ci.plac.lose.off90 <- quantile(allplac.lose.off,c(0.05,0.95),na.rm=T)

# Diff-in-disc effects:
did.off.sep # DID estimate of 0.32
ci.did.gain.off.sep # -0.33, 1.17

did.on.sep # DID estimate of 0.33
ci.did.gain.on.sep # -0.41,1.16

# Placebo effects:
plac.off # -0.23
ci.plac.gain.off 
ci.plac.gain.off90 # -1.06,0.42

plac.on # placebo est of -0.30
ci.plac.gain.on 
ci.plac.gain.on90 # -1.24,0.31



### in one specification:
fullintxDID <- lm(run_and_win_next ~ voteshare_adj + (voteshare_adj>=0) + voteshare_adj:(voteshare_adj>=0) + factor(fips) + factor(fips):(voteshare_adj>=0) + concurrent_next:(voteshare_adj>=0),
                  data=data6[which(data6$voteshare_adj>=-bw & data6$voteshare_adj<=bw),])
# stargazer(fullintxDID,omit = ".*factor\\(fips\\).*")
xtable(summary(fullintxDID)$coefficients[c(1:3,1083:1084),])



######################
##### Appendix G #####
######################
# Mechanisms

## Sample of ballots to gather for mechanisms tests:
set.seed(02139)
on_sample <- sample(data$elec_index[which(data$YearData>2000 & data$concurrent==1)],10,replace=F)
set.seed(02139)
off_sample <- sample(data$elec_index[which(data$YearData>2000 & data$concurrent==0)],10,replace=F)

# write.csv(unique(data6[which(data6$elec_index %in% c(on_sample,off_sample)),c("elec_index","fips","gov_id","city","JURIS","STATEAB","YearData","month","elecdate","concurrent")]),"ballot checking/ballotcheck_sample2.csv",row.names=F)


## Note:
## analysis of experiment data in separate file


